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Abstract 

This  thesis  is  a  parameter  study  of  the  effects  of  branching  fraction  uncertainty 
from  nuclear  decay  on  isotope  quantities  at  later  times.  Development  of  how  to  calculate 
and  accurately  draw  random  samples  of  branching  fractions  is  done.  Overall  effects  as 
well  as  isotope  quantity  uncertainty  distributions  are  also  explored.  Isotope  quantities  are 
calculated  using  exponential  moments  methods  with  transmutation  matrices.  Uncertainty 
from  both  half-lives  and  branching  fractions  is  carried  through  these  calculations  by 
Monte  Carlo  methods.  Results  of  the  study  show  that  uncertainty  from  branching 
fractions  dominates  most  isotopes  present  from  neutron  fission  after  approximately  one 
month.  Another  result  is  that  only  20%  of  isotopes  present  at  any  given  time  have 
uncertainty  from  both  half-lives  and  branching  fractions  that  are  of  the  same  order  of 
magnitude.  The  final  program  is  both  flexible  in  the  number  and  type  of  isotopes  it  can 
input  and  output,  as  well  as  computationally  efficient. 
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MODELING  RADIOACTIVE  DECAY  CHAINS  WITH  BRANCHING 
FRACTION  UNCERTAINTIES 


I.  Introduction 


Motivation 

A  method  of  calculating  radioactive  decays  with  quantifiable  uncertainty  is 
desired  to  ensure  the  utility  of  measurements.  It  is  important  to  know  quantities  of 
isotopes  after  a  given  time,  to  more  accurately  model  nuclear  fallout  or  trace  what  type 
and  amount  of  nuclear  material  was  present  originally.  Knowing  the  uncertainty  that 
goes  along  with  measurements  taken  ensures  accurate  conclusions  are  drawn.  A  key  step 
in  determining  uncertainty  is  to  identify  its  driving  sources.  The  motivation  behind  this 
work  is  to  identify  whether  one  of  these  sources  is  branching  fractions. 

Nuclear  Decay 

Decay  chains  of  radioactive  materials  depend  on  three  things,  quantity  of  the 
isotopes  initially  present,  half-lives  of  those  isotopes,  and  the  branching  fractions  of  their 
decay.  Given  initial  quantities  of  one  or  more  isotopes,  one  can  calculate  what  isotopes 
are  present  and  their  quantity  at  a  later  time.  Knowing  this  information  is  valuable  to 
understanding  the  environmental  hazards  that  result  from  radioactive  decay,  or  being  able 
to  accurately  determine  the  quantities  of  isotopes  present  at  earlier  times. 

Half-life  Uncertainty 

Half-lives  of  radioactive  materials  span  many  orders  of  magnitude.  Because  of 
this,  calculating  the  quantity  of  each  isotope  present  in  a  decay  chain  after  some  period  of 
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time  becomes  difficult.  Different  types  of  equilibrium  may  arise  due  to  long  lived 
isotopes  decaying  to  short  lived  isotopes  or  vice  versa.  Because  of  this  large  range,  half- 
life  uncertainty  can  have  a  major  effect  on  decay  chains. 

Branching  Fraction  Uncertainty 

If  an  isotope  has  more  than  one  type  of  decay  that  results  in  multiple  daughters, 
then  the  fraction  of  decays  that  result  in  each  daughter  are  taken  to  be  that  isotope’s 
branching  fractions.  The  resulting  daughters  may  or  may  not  be  radioactive,  and  their 
branching  fractions  can  vary  by  orders  of  magnitudes.  Each  of  these  fractions  has  an 
uncertainty,  and  again  because  of  the  same  reasons  laid  out  for  half-lives,  this  can  have 
an  effect  on  decay  chains. 

Decay  Example 

To  illustrate  the  effects  of  branching  fraction  uncertainty,  consider  Figure  1.  It 
shows  a  notional  diagram  of  two  initial  isotopes  (A  and  B)  and  their  decay  chains.  Under 
the  name  of  each  isotope  is  the  time  scale  of  its  half-life.  The  only  branching  in  these 
decay  chains  is  in  the  initial  decays,  which  is  not  always  the  case,  but  is  illustrated  here 
for  simplicity.  The  percentages  on  the  arrows  for  the  decay  represent  the  average 
branching  fraction.  Lastly  it  is  assumed  there  are  100  moles  of  isotope  A  initially  present 
and  10  moles  of  isotope  B. 
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Figure  1.  Diagram  of  notional  decay  chains  from  two  isotopes. 


Notice  that  all  chains  lead  to  either  stable  isotopes  G  or  J,  and  all  of  B  ultimately 
lead  to  J.  With  no  branching  fraction  uncertainty,  eventually  all  10  moles  from  B  and  40 
moles  from  A  would  become  J.  The  remaining  60  moles  of  A  would  become  G. 
Regardless  of  half-life  uncertainty,  this  would  be  the  end  result.  During  the  time  it  takes 
for  all  these  decays  to  occur,  half-life  uncertainty  would  have  an  effect  on  which  isotopes 
were  present  at  a  given  time.  Now  consider  branching  fraction  uncertainty,  and  assume 
the  40%  that  goes  from  A  — >  D  can  vary  by  ±5%.  If  45  moles  decayed  to  D  instead  of  40 
the  end  result  would  be  55  moles  of  J,  and  55  moles  of  G.  This  is  a  change  of  10%  in  the 
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amount  of  J  compared  to  when  there  was  no  branching  fraction  uncertainty.  It  not  only 
has  an  effect  on  the  end  products,  but  it  affects  when  an  isotope  may  be  present  just  like 
half-life  uncertainty.  Now  consider  the  branch  from  B  — >  D.  If  two  moles  went  to  D, 
and  only  one  went  to  E  it  changes  how  much  and  at  what  times  H  is  present.  This  is  due 
to  the  large  difference  in  half-lives  between  D  and  E.  Numerous  scenarios  can  be 
explored  just  with  this  simple  diagram.  Try  to  imagine  the  complexity  that  results  from 
having  around  800  initial  isotopes  with  thousands  of  decay  chains.  The  potential  effect 
of  branching  fraction  uncertainty  becomes  clear. 

Problem 

Model  the  time-dependent  quantities  of  all  known  isotopes  given  an  initial 
amount  of  isotopes  present.  Include  the  uncertainties  of  those  quantities  due  to  the  half- 
life  and  branching  fraction  uncertainties  of  each  individual  isotope. 

Goal 

The  goal  of  this  work  is  to  understand  how  branching  fraction  uncertainties 
contribute  to  the  overall  uncertainty  in  the  quantity  of  isotopes  present  at  a  given  time. 
This  is  done  by  comparing  or  combining  branching  fraction  uncertainties  with  the 
uncertainty  from  half-lives  to  get  a  sense  of  if  and  when  it  becomes  a  major  source  in 
quantity  uncertainty.  Quantification  of  these  uncertainties  will  help  interested 
organizations  determine  if  experiments  designed  to  measure  branching  fraction 
uncertainty  are  warranted  and  set  priorities  for  said  experiments. 
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Scope 


Because  few  of  the  branching  fraction  uncertainties  are  known,  the  scope  of  the 
project  was  to  develop  a  tool  to  complete  a  parameter  study.  Incorporating  existing 
branching  fraction  uncertainties  was  initially  explored  but  was  determined  to  be  beyond 
the  scope  of  this  study.  This  was  due  to  issues  in  weighting  and  incorporating  these 
datasets  and  the  fact  that  their  distributions  were  probably  quoted  as  normal  when  they 
are  not.  The  program  can  read  in  isotope  data  including  metastable  states,  half-life, 
branching  fraction,  and  decay  information.  It  will  also  take  in  any  number  of  initial 
isotope  quantity  datasets.  The  output  is  quantity  and  basic  statistics  to  evaluate 
uncertainty  of  any  isotope  including  metastable  states  at  any  positive  time. 

Assumptions 

A  consistent  method  of  applying  uncertainty  is  used,  because  a  majority  of  the 
branching  fraction  uncertainties  are  unknown.  For  each  isotope  that  has  more  than  one 
type  of  decay  associated  with  it,  the  branching  fraction  with  the  smallest  mean  is  chosen. 
Its  relative  standard  deviation  (standard  deviation  divided  by  the  mean  or  a/fu)  is  assumed 
to  be  equal  to  a  percentage.  This  same  percentage  is  applied  to  the  smallest  fraction  for 
all  isotopes  in  a  given  run  of  the  program,  making  the  relative  standard  deviation  of  the 
smallest  branching  fraction  a  primary  independent  variable  in  the  parameter  study. 
Because  the  standard  deviation  of  each  branch  of  a  given  isotope  is  dependent  on  the 
others,  knowing  one  allows  you  to  calculate  the  others. 

Another  assumption  is  that  there  is  no  uncertainty  in  the  initial  amounts  of  the 
isotopes  present.  Eliminating  this  variable  allowed  the  study  to  focus  on  branching 
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fraction  uncertainty  and  how  its  effects  compare  with  the  effects  from  half-life 
uncertainty.  This  feature  could  be  incorporated  into  the  program  for  fission  fragment 
yield  uncertainties  in  the  future  because  the  same  method  used  to  draw  random  samples 
for  branching  fraction  uncertainties  applies. 
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II.  Random  Sampling  of  Branching  Fractions 


The  overall  approach  uses  Monte  Carlo  methods  to  develop  quantity 
uncertainties.  For  each  Monte  Carlo  trial,  a  random  sample  from  the  uncertainty 
distributions  for  half-lives  and  branching  fractions  is  drawn.  This  method  was  chosen 
over  trying  to  carry  these  individual  uncertainties  through  the  differential  equations  that 
define  their  decay  for  its  relative  mathematical  simplicity  and  computational  efficiency. 

The  model  used  builds  upon  the  work  of  Harr  [1],  who  focused  on  half-life 
uncertainties,  and  their  effect  on  the  uncertainty  of  gamma  activity.  His  use  of  the 
exponential  moment  methods,  method  for  sampling  from  normal  distributions  for  half- 
life  uncertainty,  and  use  of  transmutation  matrices  were  leveraged.  This  chapter 
discusses  several  distributions  that  will  lead  to  an  understanding  of  what  is  necessary  to 
randomly  sample  branching  fractions.  It  then  develops  a  method  for  sampling  these 
branching  fractions,  while  avoiding  biasing  or  catastrophic  cancellation  that  can  occur 
from  normalizing. 

Bernoulli  Distribution 

The  Bernoulli  distribution  is  a  discrete  probability  function  with  only  two  possible 
outcomes.  Discrete  meaning  there  is  only  success  or  failure  for  each  trial.  The  best 
example  of  this  is  the  flipping  of  a  coin.  The  two  possible  outcomes  are  heads  or  tails. 

The  probability  for  each  outcome  must  be  in  (0,  1),  and  the  sum  of  the  two  probabilities 
must  equal  one.  If  we  assume  the  probability  is  equal  for  both  outcomes,  then  the 
probability  for  each  trial  and  for  both  heads  and  tails  is  0.5.  If  the  coin  has  a  probability 
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of  0.7  for  heads  then  the  probability  for  tails  would  be  1  -  0.7  =  0.3.  A  Bernoulli 
distribution  only  considers  one  trial  (n  =  1). 


Binomial  Distribution 

If  more  than  one  trial  is  done  with  the  same  parameters  as  the  Bernoulli 
distribution,  then  it  is  a  binomial  distribution.  If  the  probability  for  an  outcome  is  known, 
then  for  a  given  number  of  trials  the  probability  of  getting  any  number  of  successes  or 
failures  can  be  calculated  from  the  probability  mass  function  (pmf).  The  pmf  for  a 
binomial  is  given  as: 


Binomial  pmf  = 


n! 


-px(  1  —  p)n 


x!(n-x)!r  (1) 

where  n  is  the  number  of  trials  x  is  the  number  of  successes  and  p  is  the  probability  of 
success.  Consider  the  coin  example  again  where  the  probability  of  success  was  0.5  and 
now  the  number  of  trials  is  set  at  100.  There  would  be  a  7.96%  chance  of  getting  exactly 
50  successes,  and  a  1.08%  chance  of  getting  exactly  40  successes. 


Beta  Distribution 

Now  assume  the  probability  of  success  is  unknown,  and  there  are  still  only  two 
possible  outcomes.  An  experiment  is  done  with  n  trials  to  determine  the  probability  of 
success.  It  would  make  sense  that  the  more  trials  that  are  done  the  more  confidence  there 
would  be  in  this  probability.  This  is  known  as  a  conjugate  prior  distribution,  or  taking  the 
information  that  is  known  and  building  a  distribution  from  it.  The  beta  distribution  is  the 
conjugate  prior  of  the  binomial  distribution,  and  whose  probability  density  function  (pdf) 
is: 
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(2) 


Beta(x;  a,/?) 


xa  1(1— x)^  1  xa  1(l  —  x)13  1 
fit'-'Cl-ty-'dt  B(a,/?) 


where  a  and  (>  are  the  concentration  parameters  of  the  distribution,  0  <  x  <  1,  and  the 
denominator  is  the  beta  function  B(a,  /?).  As  an  example  if  no  prior  knowledge  of  the 
probability  of  success  were  known,  and  100  trials  were  performed  with  50  successes  as 
the  result.  The  beta  pdf  would  have  parameters  a  =  5 1  and  =  51,  where  a  is  equal  to 
one  plus  the  number  of  successes  and  /?  is  equal  to  the  one  plus  the  number  of  failures. 
The  resulting  pdf  is  shown  in  Figure  2  along  with  another  beta  pdf  with  parameters  a  = 
101  and  =  101.  While  they  make  look  Gaussian,  they  are  not.  What  does  make  sense  is 
that  the  highest  probability  of  success  is  right  at  0.5,  because  1/2  of  the  trials  performed 
turned  out  successful.  What  should  also  be  noted  is  that  when  the  number  of  trials  was 
doubled  to  200  the  standard  deviation  decreased,  and  the  probability  at  0.5  increased. 

This  too  makes  sense  because  you  have  twice  as  much  information  as  in  the  first 
experiment,  that  supports  that  the  probability  is  most  likely  0.5. 
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Figure  2.  pdf  of  beta  distributions  showing  difference  in  number  of  trials.  The  variance  of  the  red 
curve  is  smaller  than  the  blue  because  more  trials  were  done  and  it  has  the  same  mode. 

Now  consider  two  other  scenarios,  both  that  have  100  trials  but  one  has  65 
successes  and  the  other  has  only  5  successes.  The  resulting  pdf  distributions  are  shown  in 
Figure  3.  The  pdf  with  65  successes  still  looks  rather  Gaussian,  but  it  has  shifted  to  the 
right  to  show  a  higher  probability  of  success  around  0.6.  The  second  distribution  is  very 
close  to  zero,  and  if  is  examined  closely  the  right  tail  is  larger  than  the  left  indicating  it  is 
not  symmetric  like  a  normal  distribution. 
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Figure  3.  pdf  of  beta  distributions  with  100  trials  but  different  number  of  successes.  The  difference 
in  the  parameters  shows  how  the  mode  and  the  max  value  of  the  pdf  are  affected. 


Branching  fractions  for  each  isotope  without  uncertainty  are  documented  in  the 
National  Nuclear  Data  Center’s  NuDat  2.6  database  [2].  It  is  assumed  that  these 
branching  fractions  are  the  means  of  the  probabilities.  However,  to  solve  for  both  the  a 
and  [i  parameters,  the  standard  deviation  or  variance  of  one  of  the  branches  must  be 
known.  The  mean  and  standard  deviation  of  a  beta  distribution  are: 


a 


f^Beta 


a+p’ 


®Beta 


a(3 


(a  +  /?)2(a  +  /?  +  !) 


(3) 
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Of  note  is  that,  unlike  a  Gaussian,  the  mean  of  a  beta  is  not  the  mode.  The  simplicity  of 
these  equations  means  its  parameters  can  quickly  be  calculated.  If  the  means  of  both 
branches  and  variance  or  standard  deviation  of  one  of  the  branches  is  known,  then  its  a 
and  [)  can  be  calculated.  The  beta  distribution  can  then  be  sampled  from  and  the  other 
branch  is  determined  by  taking  one  minus  the  sample  drawn.  This  will  not  work  for 
instances  where  the  number  of  categories  is  greater  than  two,  which  require  a  different 
distribution. 

These  three  types  of  distributions  (Bernoulli,  Binomial,  and  Beta)  are  based  on 
only  having  two  possible  outcomes.  Some  isotopes  can  decay  in  more  than  two  ways,  so 
the  generalization  to  these  distributions  that  handle  multiple  possible  outcomes  is 
presented  in  the  next  sections. 

Categorical  Distribution 

The  generalization  of  the  Bernoulli  distribution  is  the  categorical  distribution.  A 
categorical  distribution  has  j  possible  outcomes.  The  probability  of  each  of  those 
outcomes  is  known  and  they  sum  to  one.  As  an  example,  isotope  A  decays  to  isotope  B 
60%  of  the  time,  isotope  C  30%  of  the  time,  and  D  10%  of  the  time.  This  means  that 
isotope  B  has  a  probability  of  0.6,  isotope  C  of  0.3,  and  isotope  D  0.1  in  the  categorical 
distribution  and  they  sum  to  exactly  1. 

Multinomial  Distribution 

The  generalization  of  the  binomial  distribution  is  the  multinomial  distribution.  As 
with  the  binomial,  more  than  one  trial  is  done  and  the  probability  of  success  for  each 
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outcome  is  known.  The  pmf  with  the  probabilities  of  getting  a  certain  number  of 
successes  is: 

L  xi 

pmf  (xlt  x2, ...  xj)  =  n\  ]”[  (4) 

i= i  l' 

where  n  —  xi  's  the  number  of  trials,  j  is  the  number  of  categories,  X;  are  the  number 
of  successes  for  each  category  (which  sum  to  n),  and  p ,  is  the  probability  of  success  for 
each  category  (which  sum  to  1).  Equation  (4)  is  equal  to  the  binomial  pmf  in  Equation 
(1),  when  j  =  2.  This  relates  to  radioactive  decay.  Assume  there  is  no  uncertainty  in  the 
branching  fractions  of  isotope  A  from  the  previous  example.  The  probability  of  seeing  a 
certain  number  of  isotope  B,  C,  or  D  from  a  known  number  of  decays  of  A  could  be 
calculated  from  the  pmf.  The  problem  is  that  the  probabilities  can  never  be  known 
exactly;  therefore,  the  conjugate  prior  to  this  distribution  is  needed.  A  conjugate  prior 
distribution  creates  a  pdf  based  on  known  data. 

Dirichlet  Distribution 

The  Dirichlet  distribution  is  the  conjugate  prior  to  the  multinomial  distribution 
and  generalization  of  the  beta  distribution.  The  Dirichlet  has  j  possible  outcomes,  takes 
known  information  and  gives  a  pdf  for  each  of  those  outcomes.  The  Dirichlet  has  j 
concentration  parameters.  Like  the  beta  distribution  each  parameter  is  equal  to  one  plus 
the  number  of  successes  witnessed.  The  pdf  of  the  Dirichlet  distribution  is  given  in 
Equation  (5). 
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(5) 


r  ^(^o)  ]  \  a~  1 

Dir(xlf  x2, ...  Xj)  alt  a2> ...  ocj )  =  ■  \\xi 

n{=1(r(a*))  \=t 

j  j 

^  oci  =  a0  and  ^  xt  =  1 

i=l 


t=l 


where  a,  are  the  concentration  parameters,  x,-  are  the  branching  fractions  in  (0,1),  and  Y  is 
the  gamma  function.  It  can  be  shown  that  the  Dirichlet  is  the  generalization  of  the  beta 
pdf  in  Equation  (2)  by  setting  j  =  2  in  Equation  (5)  and  using  the  property 
B(x,  y)  =  F (x)  F (y)  / T  (x  +  y).  Kotz  [3]  also  proves  that  the  marginal  distributions  of  a 
Dirichlet  are  beta  distributions.  A  marginal  distribution  is  a  distribution  that  defines  the 
probability  of  a  single  or  subset  of  the  possible  outcome(s).  The  beta  distributions  that 
represent  the  marginal  distributions  of  a  Dirichlet  are: 

Beta(Xj;  ait  a0  —  aj 

In  order  to  sample  each  of  these  beta  distributions,  all  concentration  parameters  of  the 
Dirichlet  must  be  known.  Similar  to  a  beta,  the  Dirichlet  parameters  can  be  calculated 
from  the  mean  and  standard  deviation  of  each  parameter,  which  is  given  in 
Equations  (7),  (8),  and  (9). 


(6) 


ai 

Mean  of  x,-  =  Ui  =  — 
a0 


Mode  ofxj  = 


gf  -  1 
«o  ~j 


ai(a0  —  a  i) 

Standard  Deviation  ofx;  =  07  =  ^ - 

ao(ao  +  1) 


(7) 

(8) 

(9) 


Just  like  a  beta,  the  means  are  assumed  to  be  the  quoted  branching  fractions.  A  closer 
look  shows  the  mean  is  a  ratio  of  the  individual  concentration  parameter  to  the  sum  of  the 
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concentration  parameters.  In  the  case  of  interest,  this  is  simply  the  branching  fraction 
mean  times  a  scalar.  It  also  shows  that  the  mode  differs  from  the  mean.  Equation  (9) 
reveals  that  as  the  sum  of  the  concentration  parameters  increases  the  standard  deviation 
decreases  approximately  inversely  with  o.q.  This  means  the  size  of  the  scalar  (sum  of  the 
concentration  parameters,  «o)  determines  the  standard  deviation  of  all  the  possibilities.  If 
all  the  means  are  known  and  the  variance  or  standard  deviation  of  one  of  the  branches  is 
known,  then  the  sum  of  concentration  parameters  can  be  calculated.  This  then 
determines  the  standard  deviation  of  the  other  branches. 

Now  expand  on  the  example  used  in  the  categorical  distribution  section.  An 
experiment  was  done  and  it  was  measured  that  out  of  100  decays  60  times  it  went  to  B, 

30  times  it  went  to  C,  and  10  times  it  went  to  D.  The  Dirichlet  distribution  would  have 
concentration  parameters  of  61,  31,  and  11.  Using  Equation  (9)  the  standard  deviation  of 
the  branch  going  to  isotope  B  would  be  0.0482.  If  the  same  experiment  was  done  but  this 
time  500  decays  were  measured  and  300  times  it  went  to  isotope  B  then  the  standard 
deviation  would  be  0.0218.  This  makes  sense  because  the  counting  statistics  are  better  in 
the  experiment  with  500  measurements  making  the  variance  smaller. 

This  can  be  seen  by  looking  at  a  3d  plot  of  our  100  measurement  example 
Dirichlet  distribution  in  Figure  4.  The  plot  is  the  pdf  which  can  only  show  two 
parameters  at  a  time.  Figure  4  shows  the  parameters  31  and  11  representing  isotopes  C 
and  D  respectively. 
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Figure  4.  pdf  of  a  Dirichlet  distribution  with  concentration  parameters  (61,  31, 11),  of  which  the 
second  two  parameters  are  shown.  The  peak  is  centered  around  x2  =  0.3  and  jc3  =  0.1,  which 
corresponds  to  the  mode  of  the  second  two  parameters  (31, 11). 


Figure  4  shows  a  single  peak  that  is  centered  near  0.3  in  the  X2  direction,  and  0.1 
in  the  X3  direction.  This  corresponds  to  the  means  of  the  parameters  31  and  1 1  in  the 
Dirichlet  distribution  (remember  the  mean  is  not  where  the  peak  is).  The  white  line  that 
cuts  through  the  middle  of  the  graph  is  a  limit  of  possible  value.  That  is,  behind  it  those 
values  are  not  possible  in  the  distribution  because  the  sum  of  the  values  must  be  equal  to 
one.  If  the  Dirichlet  concentration  parameters  are  changed  to  the  experiment  with  500 
measurements  the  3d  plot  would  look  like  Figure  5. 
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Figure  5.  3d  plot  of  a  Dirichlet  pdf,  with  concentration  parameters  (301, 151,  51),  of  which  the 
second  two  parameters  are  shown.  The  peak  is  centered  at  the  same  place  as  in  Figure  4  because 
their  modes  are  the  same,  but  is  much  narrower  because  the  sum  of  the  parameters  is  greater. 

The  peak  in  Figure  5  is  almost  in  the  same  exact  location  as  Figure  4,  but  it  is 
much  narrower  than  Figure  4’s.  This  is  expected  since  the  means  of  both  distributions 
are  the  same,  but  the  sum  of  the  concentration  parameters  in  Figure  5  is  greater  meaning 
the  variances  of  each  parameter  are  smaller,  resulting  in  a  narrower  peak. 

The  beta  and  Dirichlet  distributions  look  to  be  the  best  choice  for  randomly 
sampling  branching  fractions.  This  is  because  the  distributions  must  be  continuous  on  the 
interval  of  [0,  1],  and  the  probability  be  0  at  both  0  and  1.  If  the  probability  is  not  0  at 
both  0  and  1,  then  an  isotope  could  potentially  eliminate  itself  or  another  isotope 
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(catastrophic  cancellation).  This  eliminates  the  possibility  of  it  being  Gaussian,  because 
it  cannot  have  a  probability  equal  to  zero.  A  truncated  Gaussian  is  not  a  realistic 
possibility  either.  This  is  because  in  order  for  the  probability  to  be  0  at  0  and  1  the 
truncation  would  have  to  be  at  0  +  s  and  1  -  s  (s  -  very  small  number  or  machine 
precision).  This  leaves  a  step  value  right  at  the  boundaries  that  is  not  realistic.  The  beta 
and  Dirichlet  distributions  are  the  most  logical  choices  that  satisfy  the  above 
requirements. 

Sampling  Beta  and  Dirichlet  Distributions 

Having  chosen  the  Dirichlet  and  beta  distributions,  an  efficient  method  for 
drawing  samples  from  these  distributions  is  needed.  This  will  start  with  determining  the 
Dirichlet  concentration  parameters.  This  method  will  also  work  with  isotopes  that  have 
only  two  branches  or  j  =  2,  since  a  Dirichlet  is  the  generalization  of  a  beta  distribution. 

Consider  a  Dirichlet  distribution  where  n  \ ,  H2,  •  •  -Mj  arc  the  means  of  the  branching 
fractions,  and  ao  is  the  sum  of  the  Dirichlet  concentration  parameters.  From  Equation  (7) 
the  individual  concentration  parameters  (a,-)  become: 

Dirfaao,  n2a0, ...  Hjcc0)  (1Q) 

where  cq  =  /r ia0 ,  0  <  /q  <  1,  £(=1/q  —  1.  and  a0  >  1.  The  standard  deviation  of 
each  variable  jq,  can  also  be  calculated  from  Equation  (9)  to  be: 

(11) 

Now  let  s  be  the  index  of  the  branching  fraction  with  the  smallest  mean,  and  (p  be  the 
relative  standard  deviation  of  that  branch.  Relative  standard  deviation,  also  called 
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coefficient  of  variation,  is  defined  as  the  standard  deviation  divided  by  the  mean  (a/ju). 
The  smallest  branch  is  chosen  because  it  has  the  largest  relative  standard  deviation 
making  it  the  most  sensitive  parameter.  Using  this  we  can  solve  for  ao  in  terms  of  (p  and 


Ms- 


iMs(l-Ms) 
a0  +  l 


<PMs 


(12) 


Solve  for  ao: 


i  2  — - 1 

1  -  Ms  -  <P  Us  Ms 

«n  =  - n -  =  - n - 1 


<M2Ms 


r 


(13) 


once  ao  is  known  from  Equation  (13)  each  concentration  parameter  is  calculated  by 
multiplying  its  mean  by  ao.  Once  the  concentration  parameters  are  known  they  are 
checked  to  make  sure  none  are  greater  than  1015.  If  a  parameter  is  greater  than  1015  and 
its  mean  is  close  to  one,  the  double  precision  random  number  will  round  to  one  because  it 
does  not  have  enough  digits.  When  a  parameter  is  greater  than  1015  its  sample  is  simply 
its  mean.  In  the  current  isotope  input  data  no  concentration  parameters  meet  this 
criterion,  but  this  was  put  in  to  prevent  potential  future  errors.  This  issue  can  be  put  off 
by  using  greater  precision  at  the  cost  of  memory.  Now  each  branch  can  be  sampled  from 
its  marginal  Bctaki',;  a,,  a0  -  a,). 

To  maximize  sampling  efficiency,  several  beta  distribution  sampling  and  rejection 
schemes  were  explored.  The  first  and  easiest  to  implement  was  to  use  a  geometric 
rejection  scheme.  The  beta  pdf  is  standardized  to  unity  at  its  peak,  which  is  developed  in 
Appendix  B,  and  referred  to  as  Betas(x;  a ,  b ).  This  means  both  x  and  Betas  must  be  in 
(0,1).  A  random  number  (Ci)  is  drawn  and  Betas(Ci;  a,  b)  is  calculated.  Then  a  second 
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random  number  (&)  is  drawn  and  if  C2  <  Betas(Ci;  a,  b),  then  dj  is  the  used  as  the  random 
sample,  otherwise  the  process  is  repeated  with  two  new  random  numbers.  While  simple, 
this  method  is  not  very  efficient  and  parsimonious,  especially  if  either  beta  distribution 
parameter  is  much  less  than  the  other. 

The  other  scheme  explored  is  called  “Rejection  using  Random-Term  Sampling” 
and  is  laid  out  by  Mathews  [4].  It  uses  a  rejection  sampling  function  s(x)  that  is  a  sum  of 
positive  terms.  The  function  is  determined  by  the  user,  and  can  be  used  as  long  as  s(x)  > 
Betas(T;  a,  b )  when  0  <  x  <  1 .  Betas  is  used  again  here  because  it  will  be  shown  later  in 
this  section  to  be  more  efficient  to  calculate  than  the  actual  beta  pdf.  The  goal  is  to  pick 
terms  that  almost  equal  the  Betas,  and  that  are  easy  to  integrate  and  invert  so  that  a 
random  sample  can  be  drawn  from  the  inverted  CDF  of  s(x).  Some  examples  of 
functions  that  have  easily  inverted  CDFs  are  constants,  linears,  polynomials,  and 
hyperbolic  trigonometric  functions  (e.g.  sinh,  cosh).  For  the  cases  explored  s(x)  was 
made  a  piecewise  function  by  using  heavyside  functions  H(x;  c)  defined  as: 


HO;  c)  =  { 


0,  x  <  c 


Ll,  x  >  c  (14) 

where  c  is  the  x  value  where  the  step  function  goes  from  0  to  1.  If  n  terms  were  used  to 
make  the  sampling  function,  then  s( x)  =  s  1  (x)  +  S2O)  +  . . .  sn(x).  Making  s(x)  into  a 
piecewise  function  with  the  range  of  each  term  being  (xm,  x,)  with  heavyside  functions 
would  look  like: 


f  L 

■O)  =  ^Si(x’)H(x,Xi-1')H(x,xO 

i= 1 


(15) 


where  for  a  beta  distribution  xq  =  0,  and  xn  =  1. 
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The  first  set  of  terms  explored  was  a  combination  of  sinh  functions  on  either  side 
of  the  Betas  peak,  and  a  constant  around  the  peak.  The  issue  was  finding  a  consistent 
way  to  calculate  the  sinh  coefficients  so  it  would  be  as  close  to  Betas  as  possible.  This 
was  especially  difficult  when  either  beta  parameter  was  much  less  than  the  other, 
meaning  the  peak  was  very  close  to  either  one  or  zero.  Another  method  explored  used  a 
system  of  two  or  four  linears,  and  a  constant.  Two  linears,  one  on  each  side  of  the  peak, 
and  a  constant  around  the  peak  was  implemented  first.  Each  linear  went  from  either  zero 
or  one  to  the  inflection  point  of  Betas  on  its  respective  side  of  the  peak.  The  inflection 
points  of  Betas  were  found  by  setting  its  second  derivative  equal  to  zero  and  solving  for 
x.  Equation  (16)  shows  this  solution  and  its  derivation  can  be  found  in  Appendix  C. 

a2  +  a(b  —  1)  ±  Jab  (a  +  b  —  1) 

*  =  a2  +  (b  —  Tjb  +  a(2b  -  1)  (16) 

where  a  =  a  -  l,  b  =  (3 -  l,  and  a  and  fi  are  the  beta  distribution  parameters. 

Once  the  inflection  points  are  known,  the  area  under  each  of  the  lines  and 
constant  was  calculated,  and  then  normalized  so  their  sum  was  equal  to  one.  This  gives  a 
probability  of  selecting  each  one  of  the  functions.  A  random  number  (fi)  is  drawn  and 
determines  which  of  the  functions  will  be  sampled.  If  one  of  the  lines  is  selected 
Equation  (17)  is  used  to  sample  from  its  inverted  CDF.  The  derivation  of  Equation  (17) 
is  shown  in  Appendix  D,  and  includes  some  simplifications  when  either  so  =  0,  or  si  =  1. 

^  _ _ (s0  +  s{)%2 _ 

So  +  V^oC1  -£2)  +S1^2 

X  =  (1  —  lt)x0  +  u(x±) 
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where  s(x)  is  the  linear  sampling  function,  xo  the  lower  bound  used,  x\  the  upper  bound 
used,  so  =  s(xo),  si  =  .vU'ik  and  C2  is  the  second  random  number  drawn,  u  is  the  inverted 
value  of  the  sampling  function  CDF  with  a  range  of  [0,1].  Once  u  is  known  the  value  of 
v  can  be  calculated.  This  makes  sense  because  when  u  =  1,  x  =  *1,  and  when  u  =  0,  x  =  xo, 
which  are  the  upper  and  lower  bounds  of  the  sampling  function  chosen.  If  the  portion  of 
s(x)  randomly  chosen  is  the  constant  then  there  is  no  need  to  find  the  inverse  of  the  CDF, 
it  is  simply  u  =  C2,  and  x  is  calculated  from  the  same  equation  as  the  linear.  Now  that  x 
has  been  determined,  BctasOr;  a,  b)  and  s(x)  are  calculated.  A  third  random  number  (C3) 
is  drawn  and  if  £3  s(x)  <  BetasOc;  a,  b )  then  x  is  accepted  as  the  random  sample,  otherwise 
it  is  rejected  and  the  process  begins  again.  Shown  below  is  a  sample  algorithm  that 
generalizes  this  process  with  n  piecewise  terms. 


Function  Sample_Beta_Linear(alpha,  beta,  n)  Result(SampleX) 

Use  Kinds,  Only:  dp 
Use  Random_Numbers 
Implicit  None 

Real(dp),  Intent(In)::  alpha,  beta 

Integer,  Intent(In)::  n  Inumber  of  linear  terms  to  use 
Real(dp)::  SampleX,  y,  PDF 
Real(dp)::  a,  b,  u 

Real(dp)::  x(0:n),  s(0:n),  Area(0:n),  xi(l:3) 

Integer: :  i 

a  =  alpha  -  l._dp 
b  =  beta  -  l._dp 

IFind  x  values  where  piecewise  terms  meet 
x(0)  =  0._dp 

ForAll  (i  =  l:n-l)  x(i)  =  Calc_x_Piecewise_Points(i,  a,  b) 
x(n)  =  l._dp 
Area  =  0._dp 
Do  i  =  0,  n 

IFind  y  values  where  piecewise  terms  meet 
s(i)  =  Sampling_Function_Value(i,  x(i)) 

ICalculate  area  under  line  for  probability  of  selecting  it  to  sample  from 
If  (i  >  0)  Area(i)  =  Get_Area_Under_Curve(x(i-l) ,  x(i),  s(i-l),  s(i)) 

End  Do 

Area(0)  =  Sum(Area(l : n) ) 

ForAll  (i  =  l:n)  Area(i)  =  Area(i)/Area(0) 
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Do 


!Get  random  numbers 
xi  =  Random_Vector(3) 

Do  i  =  1,  n 

ICheck  to  see  which  line  to  sample  from  first  random  number 
If(xi(l)  <=  Sum(Area(l:i)))  Then 

! Calculate  u  from  lines  inverse  CDF  using  second  random  number 
u  =  Inverse_CDF_line(xi(2),  x(i-l),  x(i),  s(i-l),  s(i)) 

! Calculate  x  from  u,  and  y  of  sampling  function  from  x 
SampleX  =  (l._dp-u)*x(i-l)+u*x(i) 
y  =  Sampling_Function_Value(i,  SampleX) 

Exit 
End  If 
End  Do 

ICalculate  value  of  standardized  beta  PDF  from  x  drawn 
PDF  =  Beta_PDF (SampleX,  a,  b) 

ICheck  for  acceptance  by  comparing  PDF  value  to  third  random  number 
Itimes  the  value  of  the  sampling  function 
If  (xi(3)*y  <=  PDF)  Exit 
End  Do 


Using  two  linears  that  go  from  either  zero  or  one  to  the  inflection  points  and  a 
constant  does  improve  the  parsimony  of  random  sampling  over  having  no  rejection 
scheme,  however  this  improvement  is  diminished  greatly  as  either  beta  parameter 
becomes  much  greater  than  the  other.  Figure  6  shows  BctasCt;  50,  450)  and  the  sampling 
functions  with  three  terms  as  an  example  of  when  this  happens. 
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Figure  6.  Betas(x;  50,  450),  orange  and  red  lines  are  rejection  sampling  functions  from  0  to  the  1st 
inflection  point,  and  2nd  inflection  point  to  1  respectively.  Green  line  represents  constant  sampling 

function  around  the  peak  of  Betas. 

What  is  apparent  from  Figure  6  is  that  if  the  areas  under  the  three  sampling 
functions  are  compared,  the  area  under  the  red  curve  dominates.  This  means  the  red  line 
will  be  chosen  to  sample  from  a  majority  of  the  time.  The  problem  is  the  value  of  Betas 
under  the  red  line  is  almost  zero  making  the  probability  of  the  third  random  number  being 
less  than  the  Betas  very  small.  This  highlights  the  need  for  another  set  of  lines  that  start 
at  0  and  1  and  go  to  points  very  close  to  where  the  peak  begins  to  rise.  Consistently 
being  able  to  find  these  points  over  a  wide  range  of  beta  parameters  is  the  challenge. 
Equations  (18)  and  (19)  were  empirically  determined  and  represent  these  points,  the 
derivation  of  which  is  shown  in  Appendix  E.  As  in  development  of  the  beta  distribution 


24 


inflection  point,  a  =  a  -  1  and  b  =  J3  -  1,  where  a  and  f>  are  the  beta  distribution 
parameters. 

cl  +  b  +  SR  — 

2S^  08) 

where  Sr  is  the  relative  slope  (Betas'(x;  a,  b)/ Betas(x;  a,  h))  to  the  left  and  right  of  the 
peak  are: 

c  ±(a  +  h) 

~  el°gio(a)-l  (19) 

Now  a  system  of  four  lines  and  one  constant  can  be  used.  The  first  line  goes  from  zero  to 
the  positive  relative  slope  point  in  Equation  (18).  The  second  line  goes  from  this  point  to 
the  first  inflection  point.  The  constant  covers  the  area  between  the  two  inflection  points. 
The  third  line  goes  from  the  second  inflection  point  to  the  negative  relative  slope  point. 
The  last  line  covers  from  the  negative  relative  slope  point  to  one.  Figure  7  shows  the 
same  Betas(x;  50,  450)  from  Figure  6,  but  with  this  new  set  of  sampling  functions. 
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Figure  7.  Betas(50,  450),  orange,  and  red  lines  represent  the  four  linear  rejection  scheme.  Green  line 
represents  constant  sampling  function  around  the  peak  of  Betas. 

The  two  lines  that  cover  from  zero  and  one  to  the  relative  slope  points  can  barely 
be  seen  in  Figure  7,  but  neither  can  Betas.  This  makes  the  four  linear  sampling  scheme 
much  more  desirable,  because  the  area  under  the  two  relative  slope  curves  is  very  small 
compared  to  the  constant  and  other  two  lines.  This  means  the  relative  slope  lines  are 
selected  far  less  to  sample  from,  and  when  they  are  sampled  there  is  a  higher  probability 
the  sample  will  be  accepted.  If  efficiency  is  considered  to  be  the  area  under  the  Betas 
curve  compared  to  the  area  under  s(x)  then  Figure  7  is  much  more  efficient  than  Figure  6. 
The  four  linear  and  constant  scheme  was  tested  for  a  wide  range  of  beta  parameters 
ranging  from  1  to  1015.  All  test  cases  ended  up  having  acceptance  rates  from  50%  to 
75%,  with  a  majority  being  in  the  low  60%  range. 
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Now  that  an  effective  rejection  scheme  has  been  developed,  an  efficient  method 
for  calculating  the  beta  pdf  is  needed.  Recall  the  equation  of  the  beta  pdf  with  parameters 
a  and  /?  from  Equation  (2). 


Beta(x;  a,/?) 


xa  X(1  —  x)P  1 
B  (a,/?) 


where  B(a,/?)  =  f^ta  X(1  —  t)^  1dt 

There  are  two  problems  with  this  form  of  the  equation.  One,  with  the  integral  in  the 
denominator,  it  makes  the  calculation  computationally  expensive  due  to  the  need  for  a 
table  look  up  or  quadrature  routine.  Two,  the  distribution  peak  value  varies  depending  on 
a  and  [>.  These  two  issues  are  solved  by  standardizing  the  distribution  to  where  the  peak 
value  is  equal  to  one.  This  is  done  by  dividing  Equation  (2)  by  its  peak  value.  This  gives 
a  standardized  pdf  (Betas(x;  a ,  b ))  whose  values  have  a  range  of  [0,1]  and  eliminates  the 
B(a,  [:>)  integral  through  cancellation.  Doing  this  along  with  the  simplification  of  using  a 
=  a  -  1,  and  b  =  J3  -l,  yields  Equation  (20). 

(cl  +  b^)a+b 

Betas (x;  a  =  a  -  l,b  =  p  -  1)  =  — — xa(l  -  x)b  (20) 


Equation  (20)  is  now  computationally  inexpensive  to  calculate.  A  new  issue  arises  with 
the  exponent  ( a  +  b )  in  the  numerator.  Numerical  overflow  can  now  occur,  and  testing 
revealed  it  does  using  double  precision  when  a  +  b  >  143.  This  limit  is  not  acceptable,  as 
a  and  fi  values  from  Dirichlet  distributions  for  isotopes  can  be  as  large  as  10  .To  solve 
the  overflow  issue  logarithms  are  used.  Now  Betas  can  be  represented  by  Equation  (21) 
without  concern  for  numerical  overflow. 


Betas(x;a,  b)  =  e (a+fc>) log(a+i>)— ot log(ci)— fc>  log(fo)+ci log(^)+fc>  log(l-x) 


(21) 
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Equation  (21)  uses  five  logarithms  and  one  exponential  evaluations  rather  than  five  and 
five  in  Equation  (20).  Equation  (21)  is  expected  to  take  only  60%  as  long  to  compute. 
This  was  confirmed  in  testing  over  a  range  of  a  and  /?  values  that  did  not  overflow 
Equation  (20).  A  derivation  of  all  three  equations  can  be  found  in  Appendix  B.  With  a 
well-conditioned  and  relatively  computationally  inexpensive  method  for  calculating  the 
value  of  Betas,  all  the  pieces  are  in  place  to  draw  random  samples  for  branching  fractions. 
This  allows  for  calculation  of  isotope  quantities  at  later  times  in  an  efficient  manner. 
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III.  Implementation 


The  program  is  divided  up  into  six  main  modules:  reading  in  the  problem  and 
input  data,  mapping  and  organizing  the  isotope  and  gamma  data,  sampling  from 
distributions  for  half-life  and  branching  fraction  uncertainties,  building  decay  chains  and 
generating  the  transmutation  matrix  (T-matrix),  calculating  isotope  quantities  using  the 
exponential  moment  functions,  and  finally  calculating  statistics  and  writing  the  output 
files.  This  chapter  will  lay  out  details  for  each  of  these  modules  except  for  decay 
chain/T-matrix  generation,  and  exponential  moments  calculations.  These  modules  only 
have  minor  changes  from  Harr’s  [1]  code  such  as  variable  name  changes  and  interfaces 
with  other  subroutines/modules.  A  flowchart  that  includes  how  all  the  subroutines  and 
functions  interact  and  are  organized  into  modules  can  be  found  in  Appendix  A. 

Code  Design 

The  code  was  written  in  Fortran  95  using  the  Microsoft  Visual  Studio  2010 
development  suite.  Fortran  was  chosen  for  its  speed  in  performing  mathematical 
calculations,  and  simple  programming  language.  Several  helper  modules/subroutines  are 
included  in  the  code.  These  include  sorting  routines,  random  number  generators,  timing 
routines,  a  variables  module,  a  files  module,  and  a  kinds  module.  The  sorting  routines 
and  random  number  generators  are  self-explanatory.  The  timing  routines  calculate  the 
exact  time  in  seconds  in  linear  or  log  time  (specified  in  the  input  data),  as  well  as  the  total 
number  of  times  to  be  run.  The  variables  module  holds  all  the  global  variables  of  any 
type  used  in  the  code.  This  makes  it  easy  to  find  information  on  any  variable  that  is  used 
in  multiple  routines,  and  minimized  parameters  that  need  to  be  passed  between  routines 
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and  functions.  The  files  module  holds  an  organized  list  of  named  integers  that  represent 
the  extemal/internal  unit  specifier  for  writing  and  reading  files.  This  is  helpful  to  ensure 
that  the  same  unit  specifier  is  not  inadvertently  used  in  different  routines  when  files  are 
still  open.  The  kinds  module  has  specifiers  for  single,  double,  and  quadruple  machine 
precision  for  numbers.  Currently  the  entire  program  uses  double  precision  and  is  able  to 
sample  from  distributions  and  perform  calculations  without  loss  of  precision.  There  may 
be  cases  in  the  future  when  a  small  enough  branching  fraction  is  measured  that  would 
require  quadruple  precision.  Another  possible  need  for  greater  precision  arises  if 
uncertainties  in  fission  fragment  yields  are  sampled. 

Input  Data 

The  first  information  read  into  the  program  is  the  problem  data.  This  file  holds 
the  names  of  the  isotope  and  gamma  data  files,  as  well  as  the  initial  isotope  formation 
files.  It  also  defines  which  uncertainties  will  be  used,  half-life  and/or  branching  fractions. 
As  a  holdover  from  the  Harr  [1]  code,  activity  calculations  and  their  outputs  is  also  an 
option.  Other  options  include  linear  or  logarithmic  time  steps,  and  their  start  and  stop 
time  or  decades.  The  value  of  <p,  discussed  in  the  Assumptions  section  is  also  set  in  this 
file.  Some  other  performance  options  are  the  number  of  Monte  Carlo  trials  to  perform, 
the  maximum  number  of  decays  per  chain,  and  the  line  number  in  the  isotope  data  of  the 
highest  parent  isotope  present.  The  maximum  number  of  decays  per  chain  has  to  do  with 
array  allocation,  thus  memory  use.  With  the  current  data  being  used  the  longest  decay 
chain  is  25  isotopes  long.  This  is  a  good  number  to  start  with,  however  the  option  is  put 
in  the  problem  data  because  as  new  isotope  data  is  published  decays  chains  can 
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potentially  grow.  The  highest  parent  isotope  present  will  have  the  program  only  build  the 
T-matrix  and  do  quantity  calculations  up  to  this  specified  isotope.  This  can  dramatically 
increase  the  speed  of  the  program,  but  is  dependent  on  the  initial  isotopes  present.  For 
the  different  types  of  neutron  induced  fission  used  in  this  work,  the  highest  isotope 
initially  present  was  Lutetium  172.  As  a  safety  margin  the  highest  parent  isotope  was 
input  as  Francium  233.  This  increased  program  efficiency  by  -15%  over  calculating  all 
known  isotopes. 

The  next  file  read  into  the  program  is  the  isotope  data  file.  This  comes  from  the 
National  Nuclear  Data  Center’s  NuDat  2.6  database  [2].  The  specific  format  used  comes 
from  the  Nuclear  Wallet  Cards  output  with  uncertainties  in  the  Nuclear  Data  Sheets  style. 
The  data  is  first  sorted  by  number  of  protons  Z,  then  by  number  of  neutrons  N,  then  by 
meta  stable  state  or  Energy,  and  lastly  by  branching  fraction  percent.  All  columns  are 
read  in,  but  only  the  Z,  N,  Energy,  Tl/2  (txt),  Tl/2  (seconds),  Dec  Mode  (decay  mode), 
and  Branching  (%)  columns  are  used  by  the  program.  A  sample  of  the  data  format  is 
shown  in  Table  1,  and  skips  from  isotope  H6  to  Na24  to  show  a  metastable  state  via  the 
Energy  column. 
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Table  1.  Sample  Nuclear  Wallet  Card  output  from  the  NuDat  2.6  database. 


A 

Element 

Z 

N 

Energy 

Mass  Exc 

Unc 

Tl/2  (txt) 

Tl/2 

(seconds) 

Abund. 

Unc 

Dec 

Mode 

Branching 

(%) 

1 

N 

0 

1 

0 

8.0713 

5 

10.183  M  17 

610.98 

-1.00E+03 

-1000 

B- 

100 

1 

H 

1 

0 

0 

7.2889 

10 

-1000 

STABLE  NA 

-1000 

1.00E+00 

70 

X 

-1000 

2 

H 

1 

1 

0 

13.136 

15 

-1000 

STABLE  NA 

-1000 

1.15E-04 

70 

X 

-1000 

3 

H 

1 

2 

0 

14.9498 

22 

12.32  Y  2 

388800000 

-1.00E+03 

-1000 

B- 

100 

5 

H 

1 

4 

0 

32.89 

9 

5.7  MEV21 

7.995E-23 

-1.00E+03 

-1000 

2N 

100 

6 

H 

1 

5 

0 

41.9 

3 

1.6  MEV  4 

2.848E-22 

-1.00E+03 

-1000 

N 

100 

24 

Na 

11 

13 

0 

-8.4179 

4 

14.997  H  12 

53989.2 

-1.00E+03 

-1000 

B- 

100 

24 

Na 

11 

13 

0.4722 

-7.9457 

4 

20.18  MS  10 

0.0202 

-1.00E+03 

-1000 

IT 

99.95 

24 

Na 

11 

13 

0.4722 

-7.9457 

4 

20.18  MS  10 

0.0202 

-1.00E+03 

-1000 

B- 

0.05 

25 

Na 

11 

14 

0 

-9.3578 

12 

59.1  S  6 

59.1 

-1.00E+03 

-1000 

B- 

100 

The  next  file(s)  to  be  read  in  is  the  initial  isotope  quantity  data.  These  files 
contain  the  Element,  A  number,  metastable  state,  and  initial  quantity.  A  sample  of  the 
format  is  shown  in  Table  2.  Multiple  files  with  different  initial  quantity  information  can 
be  used.  As  will  be  seen  later  in  the  Output  Statistics  and  Data  section  the  different  initial 
amounts  will  have  their  own  rows  or  columns  in  all  the  different  output  files. 


Table  2.  Sample  format  for  the  initial  isotope  quantity  input. 


A 

Element 

Meta 

Amount 

67 

Mn 

0 

1.58E-09 

67 

Ni 

0 

1.17E-06 

67 

Zn 

0 

1.92E-10 

68 

Co 

0 

2.23E-06 

68 

Cu 

1 

3.76E-07 

68 

Cu 

0 

1.39E-07 

68 

Fe 

0 

1.61E-07 

The  initial  isotope  quantities  used  for  this  study  were  eight  types  of  neutron 
induced  fission.  They  were  U235  fission  with  thermal,  fast,  and  high  energy  neutrons; 
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U238  fast,  and  high;  as  well  as  Pu239  thermal,  fast,  and  high  energy  neutrons.  The  data 
used  for  these  different  fissions  comes  from  England  and  Rider  [5].  Neutron  induced 
fission  was  chosen  because  it  realistically  yields  the  largest  number  of  initial  isotopes 
(average  of  821).  This  brings  into  play  as  many  decay  chains,  and  initial  quantities  at 
different  points  in  the  chains  as  possible. 

Once  the  initial  quantities  have  been  read  in,  the  isotopes  of  interest  file  is  read  in 
to  the  program.  This  file  lists  all  the  isotopes  for  which  output  files  will  be  written.  This 
could  be  from  one  to  all  known  isotopes.  The  reason  for  this  file  is  if  all  isotopes  are 
output  it  will  take  more  time  to  write  all  the  output  files,  and  will  most  likely  take  up  in 
the  gigabyte  range  of  hard  drive  space.  If  the  user  is  interested  in  all  isotopes,  a  separate 
program  was  written  that  can  read  in  the  same  isotope  data  file  and  output  an  isotope  of 
interest  file  that  includes  all  isotopes.  It  is  recommended  that  if  the  isotope  data  file  is 
modified  or  a  new  one  is  created,  to  run  it  through  this  other  program  to  ensure  all 
potential  isotopes  are  included. 

The  last  files  to  be  read  are  those  needed  if  activity  is  going  to  be  calculated. 
These  files  for  gamma  data,  and  gamma  bin  information  are  in  the  format  used  by  Harr 
[1].  Asa  final  note  for  this  section,  all  the  input  files  except  for  the  problem  data  file 
require  an  extra  header  line  that  gives  either  the  total  number  of  lines  in  the  file,  or 
number  of  isotopes  of  interest,  number  of  initial  quantity  files,  etc.  After  any  changes  to 
a  file  this  number  needs  to  be  checked  before  the  program  is  run. 


33 


Isotope  Data 

Once  all  the  input  data  has  been  read  in,  some  checking,  manipulation,  and 
mapping  of  the  isotope  data  needs  to  be  done.  The  first  is  a  branching  check  of  all 
isotopes.  Before  a  check  is  done  to  make  sure  each  isotope’s  different  decays  add  up  to 
one,  some  things  are  changed.  For  instance  in  cases  where  an  isotope  undergoes  B- 
decay  and  something  else  at  the  same  time  the  NuDat  data  includes  that  other  portion  of 
the  decay  in  the  B-  branching  fraction.  As  an  example  let’s  say  an  isotope  undergoes  just 
B-  decay  95%  of  the  time,  and  B-  plus  alpha  decay  5%  of  the  time.  The  NuDat  data  will 
show  two  lines  one  with  B-  as  100,  and  BA  (for  B-  plus  alpha)  as  5,  the  fractions  of 
which  do  not  sum  to  one.  The  branching  fraction  check  subroutine  corrects  issues  like 
these  by  simply  looking  for  the  subsets  of  B-  and  electron  capture  decay  and  subtracting 
them  accordingly.  After  this  is  done  the  decay  fractions  of  each  isotope  is  summed  and  if 
it  is  not  within  a  set  tolerance  of  one,  the  line  number,  Z,  N  and  branching  fraction  sum  of 
the  isotope  is  written  to  an  error  file. 

After  the  branching  fraction  check,  if  half-life  uncertainty  is  used  a  subroutine 
will  take  the  text  portion  of  the  isotope  data  input  file  and  extract  the  standard  deviation 
and  turn  it  into  a  real  number.  There  are  cases  where  the  standard  deviation  is  quoted  as 
different  numbers  above  and  below  the  mean.  In  this  case  it  will  calculate  two  standard 
deviations,  one  that  is  used  for  values  above  the  mean  and  one  for  below  the  mean. 

The  last  thing  done  to  the  isotope  data  is  to  map  each  decay  to  its  daughter 
isotope.  This  is  done  by  going  through  each  decay  from  the  isotope  input  data  file  and 
selecting  the  case  it  matches.  That  case  will  determine  how  many  protons,  neutrons,  or 
metastable  states  to  add  or  subtract  from  the  parent  isotope.  It  will  then  find  where  this 
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daughter  isotope  is  within  the  isotopes  array  and  put  a  pointer  for  that  decay  to  it.  If  a 
case  for  the  decay  is  not  found  it  is  written  to  an  error  output  file.  Currently  a  majority  of 
the  decay  type  errors  are  spontaneous  fission.  The  program  does  not  currently  deal  with 
spontaneous  fission.  This  is  not  an  issue  because  all  of  the  initial  isotopes  formed  for  this 
study  and  their  daughters  do  not  have  enough  protons  and  neutrons  to  fission 
spontaneously. 

Once  these  routines  are  done  with  the  isotope  data,  similar  routines  go  through  the 
gamma  data  per  Harr  [1],  Once  this  is  done,  initialization  of  the  program  is  complete. 

The  first  Monte  Carlo  trial  is  started  and  uncertainties  are  sampled  from  their  respective 
distributions. 


Distribution  Sampling 

If  half-life  uncertainty  is  a  selected  option,  then  a  sample  for  each  isotope  is 
drawn  from  a  normal  distribution.  This  is  done  by  drawing  12  random  numbers, 
summing  them  and  subtracting  6.  If  the  isotope  has  two  standard  deviations,  the  one  for 
below  the  mean  is  used  if  the  value  is  negative,  and  the  one  for  above  the  mean  is  used  if 
the  value  is  positive.  The  value  is  then  multiplied  by  the  standard  deviation  and  then 
added  to  the  mean  to  get  the  random  sample. 


t HL  —  t1 HL  +  aHL 


(  n  \ 

\i=i  / 


(22) 


where  tHi  is  the  randomly  sampled  half-life,  uhl  and  oHl  are  the  half-life  mean  and 
standard  deviation  respectively,  £  is  a  random  number  in  (0,  1),  and  n  is  the  number  of 
random  numbers  drawn  (12  in  this  case).  This  method  is  not  very  parsimonious  with 
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random  numbers,  but  is  efficient  computationally  and  covers  +/-  six  standard  deviations 
from  the  mean.  It  was  also  done  because  as  Harr  [1]  discovered  by  contacting  the 
National  Institute  of  Standards  and  Technology  (NIST),  isotopes  with  these  two  standard 
deviations  represented  normal  distributions  with  different  standard  deviations  on  each 
side  of  the  peak. 

If  branching  fraction  uncertainty  is  selected  then  the  method  laid  out  in  the 
Sampling  Beta  and  Dirichlet  Distributions  section  is  used.  This  method  is  not  as 
straightforward  as  the  half-life  uncertainty  sampling  method.  First,  each  isotope  is 
checked  to  see  if  it  has  multiple  branches.  If  there  are  multiple  branches  then  a  Dirichlet 
sampling  routine  is  called  with  the  isotope  information  and  number  of  branches.  Based 
on  the  means  of  each  branch  and  the  input  parameter  (p,  the  sum  of  the  concentration 
parameters  is  calculated  and  then  each  concentration  parameter  individually.  From  there 
a  subroutine  that  implements  the  rejection  scheme  developed  earlier  is  called  to  get  a 
random  sample  for  each  branching  fraction.  After  all  branches  of  an  isotope  are 
randomly  sampled,  they  are  normalized  so  their  sum  is  equal  to  one.  This  process  is  then 
repeated  for  all  isotopes. 

Output  Statistics  and  Data 

The  last  piece  of  the  program  takes  all  the  quantity  calculations  and  performs 
some  statistical  calculations  and  outputs  the  data  to  text  files  for  analysis.  After  each 
Monte  Carlo  trial  a  routine  is  called  that  writes  the  quantity  values  for  each  isotope  of 
interest  to  a  file.  When  complete  there  is  a  separate  file  for  each  isotope  of  interest  that 
contains  all  times  run,  all  initial  quantity  inputs,  and  the  quantity  values.  If  there  are 
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numerous  isotopes  of  interest,  writing  all  the  quantity  values  to  a  file  will  take  up  a  large 
amount  of  memory,  as  well  as  time.  In  the  cases  run  for  this  study  that  did  1,000  Monte 
Carlo  trials,  with  31  different  times,  and  8  initial  quantity  inputs  each  isotope  had  a  file 
size  of  approximately  4.5  Megabytes.  With  all  known  isotopes  input  as  isotopes  of 
interest  it  takes  up  about  3.6  Gigabytes  of  disk  space  for  each  run.  If  an  isotope  of 
interest  has  quantity  values  of  zero  at  all  times  run,  for  all  initial  quantity  inputs  then  a 
file  will  not  be  created  for  it.  This  is  done  as  to  not  waste  disk  space,  and  computer  time 
writing  out  zeros  to  a  file.  Because  of  the  potential  need  for  large  amounts  of  disk  space, 
writing  out  all  these  quantity  values  is  an  option  that  is  specified  in  the  problem  data  file. 

Once  all  the  Monte  Carlo  trials  are  completed,  the  mean  and  standard  deviation  of 
each  isotope  of  interest  is  calculated  from  its  quantity  values.  If  the  isotope  of  interest 
has  values  greater  than  zero  at  any  time  run  then  two  files  are  opened  one  for  its  mean 
and  the  other  for  standard  deviation.  Its  respective  mean  and  standard  deviations  are  then 
written  to  these  files  for  each  time  and  initial  quantity  input. 

Another  option  in  the  problem  data  is  to  output  the  overall  statistical  data.  If  this 
option  is  marked  true  then  the  relative  standard  deviation  of  each  isotope  of  interest 
present  is  calculated.  Then  the  average  and  standard  deviation  of  these  relative  standard 
deviations  are  calculated  for  each  time  run  and  initial  quantity  input.  The  results  are 
written  to  two  files.  A  third  file  is  written  containing  a  list  of  the  relative  standard 
deviations  with  isotopes  from  greatest  to  least  for  each  time  and  initial  quantity.  This 
allows  the  user  to  see  which  isotopes  of  interest  have  the  largest  relative  standard 
deviations. 
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Verification 


To  ensure  the  overall  code  performed  as  expected  the  modules  themselves  and 
even  routines/functions  were  tested  individually.  More  in  depth  verification  was  done 
specifically  for  beta  distribution  sampling.  A  Monte  Carlo  program  was  written  that 
would  randomly  sample  a  beta  distribution  with  and  without  rejection  schemes.  Test 
cases  had  distribution  parameters  that  ranged  from  1  to  1024,  with  any  run  that  had 
parameters  greater  than  1015  using  quadruple  precision.  Because  random  samples  were 
being  taken  the  outputs  were  not  exact.  For  test  cases  performed  with  1,000  Monte  Carlo 
trials  the  mean  and  standard  deviation  only  varied  by  one  hundredth  of  a  percent. 

Beta  pdf  values  were  verified  by  using  the  same  range  of  parameters  above  and 
comparing  results  from  Equation  (21)  to  results  from  Equation  (20).  This  proved  that 
Equation  (21)  would  not  overflow  while  returning  an  answer  that  had  the  same  precision 
as  Equation  (20)  to  13  digits. 

The  overall  isotope  quantity  outputs  were  verified  by  comparing  the  output  of  the 
program  without  half-life  or  branching  fraction  uncertainty  to  a  separate  program  that 
calculated  these  same  quantities  using  both  a  numerical  integration  solution  and  a 
Bateman  equation  solution.  The  outputs  between  the  Bateman  solution  and  the  program 

o 

(exponential  moments  solution)  typically  had  maximum  relative  errors  of  around  10'  . 
The  relative  error  between  the  numerical  integration  solution  and  the  program  had 
maximums  of  10'“. 
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Performance 


All  runs  and  development  were  done  on  a  Dell  XPS  8300  with  an  Intel  i7-2600 
processor  at  3.40  GHz,  and  12  gigabytes  of  RAM.  The  operating  system  is  Windows  7 
Professional  64  bit.  Service  Pack  1.  The  development  suite  used  was  Microsoft  Visual 
Studio  2010  version  10.0.30319.1,  with  the  Intel  Visual  Fortran  Composer  XE  2011 
Update  9,  version  12.1.3526.2010.  The  program  can  be  compiled  to  run  using  either  32 
or  64  bit.  Initially  Intel’s  auto  parallelization  mode  was  selected  when  compiling  but 
after  testing  it  was  found  disabling  it  actually  decreased  run  times  by  about  75%. 

Typical  run  times  using  this  setup  along  with  1,000  Monte  Carlo  trials,  31  time 
steps,  all  known  isotopes  as  isotopes  of  interest,  and  8  different  initial  quantity  inputs 
were  around  2.5  hours.  The  run  time  is  most  sensitive  to  the  number  of  Monte  Carlo 
trials  and  the  number  of  time  values.  Program  initialization  takes  about  45  seconds,  and 
writing  final  statistics  to  files  takes  only  a  few  seconds.  Each  Monte  Carlo  trial  has  an 
overhead  computation  time  of  about  one  second,  which  is  primarily  the  time  needed  to 
sample  all  the  half-life  and  branching  fraction  uncertainties.  This  means  each  time  step, 
or  the  time  it  takes  it  to  generate  a  transmutation  matrix  and  calculate  all  the  exponential 
moments,  takes  about  1/4  of  a  second. 

The  same  setup  uses  approximately  150  megabytes  of  RAM  while  running.  This 
relatively  low  memory  requirement  is  achieved  by  writing  out  the  individual  isotope 
quantity  amounts  after  each  Monte  Carlo  run.  If  this  was  not  done  the  same  run  would 
require  almost  5  gigabytes  of  RAM  during  a  run.  This  is  an  issue  for  two  reasons.  One 
32  bit  Windows  systems  cannot  handle  arrays  that  require  more  than  2  gigabytes  of 
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memory.  Two  it  places  a  large  strain  on  a  majority  of  PCs  currently  in  use,  and 
eliminates  the  possibility  of  doing  two  runs  on  the  same  computer  at  once. 

As  mentioned  earlier  the  time  to  sample  all  the  branching  fractions  and  half-lives 
currently  takes  about  one  second,  which  is  about  1 1%  of  the  computer  time.  This  is  with 
the  sampling  rejection  scheme  that  uses  four  linears  and  a  constant.  If  no  rejection 
scheme  is  implemented  the  number  of  random  numbers  used  explodes  and  the  sampling 
routine  dominates  the  program  runtime.  To  test  the  rejection  scheme’s  performance  a 
program  run  with  the  same  inputs  as  one  of  the  analysis  runs  was  done.  Counters  were 
added  to  record  the  number  of  calls  to  the  beta  pdf  sampling  routine  as  well  as  the 
number  of  random  numbers  used.  For  one  Monte  Carlo  trial  without  a  rejection  scheme 
the  beta  pdf  was  called  2,179  times,  and  used  5,997,231,659  random  numbers.  That  is  an 
average  of  2,752,286  random  numbers  per  call  to  the  routine.  The  run  took  53  minutes, 
which  is  over  a  3,000  times  increase  in  run  time.  With  the  rejection  scheme  the  sampling 
routine  was  called  the  same  number  of  times  (2,179),  but  used  only  10,464  random 
numbers.  That  is  an  average  of  4.8  random  numbers  for  each  call  to  the  routine.  The 
minimum  number  of  random  numbers  the  routine  uses  is  three,  which  means  it  has  an 
average  acceptance  rate  of  63%.  This  proves  the  need  for  and  effectiveness  of  the  four 
linear  and  constant  rejection  scheme  over  a  wide  range  of  beta  distribution  parameters. 

Program  run  times  are  on  the  order  of  1  to  3  hours  for  1,000  Monte  Carlo  trials. 
Further  improvements  could  be  found  in  how  the  transmutation  matrices  are  generated. 
The  next  step  that  would  most  likely  yield  the  most  improvement  in  run  time  is  to  save  all 
the  decay  chains  to  an  array,  so  it  only  has  to  be  calculated  once. 
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IV.  Results  and  Analysis 


To  investigate  the  effect  of  branching  fraction  uncertainty  on  isotope  quantities, 
four  branching  fraction  relative  standard  deviation  percentages  (defined  as  <p  in  the 
Sampling  Beta  and  Dirichlet  Distributions  section)  were  explored:  1%,  5%,  10%,  and 
20%.  For  each  of  these  (ps ,  two  separate  runs  were  completed,  one  where  only  branching 
fraction  uncertainties  were  used,  and  another  where  both  branching  fraction  uncertainties 
and  half-life  uncertainties  were  considered.  A  single  run  was  done  with  only  half-life 
uncertainties,  because  it  is  independent  of  (p.  Short  hand  notation  is  used  hereafter  for 
each  type  of  run  done,  BF  for  branching  fraction  uncertainty  only,  HL  for  half-life 
uncertainty  only,  and  HLBF  for  runs  done  with  both  these  uncertainties. 

The  fission  fragment  yields  from  three  neutron  spectrums  of  three  different 
isotopes  were  used  as  initial  quantity  values.  The  neutron  spectrums  were  thermal,  fast 
(fission  energy  neutrons),  and  high  (fusion  energy  neutrons).  The  three  fuel  isotopes 
were  U235,  U238,  and  Pu239  (note:  U238  does  not  have  a  thermal  fission  cross  section, 
so  no  data  for  this  case  is  presented).  This  makes  eight  initial  quantity  values  that  will 
henceforth  be  referred  to  by  their  fuel  isotope  then  neutron  spectrum  (e.g.  Pu239  high  or 
U235  thermal). 

The  time  values  were  logarithmic  starting  at  1 0  1  seconds  and  ending  at  109 
seconds  (-31.7  years),  with  three  time  values  per  decade.  The  range  was  chosen  because 
it  represents  times  of  interest  for  nuclear  forensics  and  fallout  prediction,  and  was  a  large 
enough  range  to  recognize  and  analyze  trends.  1,000  Monte  Carlo  trials  were  performed 
for  two  runs  with  the  same  input  parameters  to  compare  relative  error.  This  was  done  to 
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see  if  1,000  Monte  Carlo  trials  was  a  statistically  large  enough  sample.  Between  the  two 
runs,  the  mean  isotope  quantities  had  an  average  relative  error  of  around  10  4.  The 
standard  deviation  average  relative  error  was  10'".  This  shows  the  means  typically  agreed 
to  four  digits,  and  the  standard  deviations  to  two  digits.  Two  longer  runs  with  10,000 
Monte  Carlo  trials  were  performed  and  again  for  the  mean  the  average  relative  error  was 
10'  .  The  standard  deviation  relative  error  was  10  ".  Because  increasing  the  number  of 
trials  to  10,000  had  no  or  very  little  effect  on  relative  error,  1,000  trials  was  used  for  the 
remaining  runs. 

Uncertainty  Ratios 

In  order  to  compare  and  analyze  the  quantity  of  all  isotopes  present,  a  metric  that 
is  independent  of  magnitude  is  needed.  This  is  because  the  amount  of  an  isotope  can 
vary  by  orders  of  magnitude  over  a  period  of  time,  but  also  the  range  of  quantities 
between  isotopes  varies  by  orders  of  magnitude.  Relative  standard  deviation  (cr//r)  is 
properly  suited  for  this,  but  in  this  case  its  overall  statistics  get  skewed  to  larger  numbers 
from  those  isotopes  that  are  rapidly  decaying  away  to  zero  quantity.  When  this  happens 
to  an  isotope  its  relative  standard  deviation  gets  very  large,  potentially  greater  than  one. 

Because  of  the  relative  standard  deviation  problem,  ratios  of  the  different  standard 
deviations  were  used.  For  each  isotope  present  at  a  given  time  the  ratio  of  its  standard 
deviation  from  the  BF  run  was  divided  by  the  standard  deviation  from  the  HLBF  run 
(obf/chlbf)-  Similarly  the  standard  deviation  from  the  HL  run  was  divided  by  the 
standard  deviation  of  the  HLBF  run  (oHl/ oHlbf)-  This  gives  the  fraction  of  uncertainty  in 
an  isotope’s  quantity  due  to  branching  fractions  and  half-lives  respectively.  The  ratios 
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will  be  in  [0,1],  because  total  uncertainty  must  increase  or  remain  the  same  when 
additional  source  of  uncertainty  is  introduced.  It  also  avoids  the  relative  standard 
deviation  problem  because  the  mean  is  not  used.  If  the  statement  above  about  the  two 
ratios  was  accurate  then  gBF  =1  — which  is  not  the  case.  These  uncertainties  are 

°HLBF  °HLBF 

not  additive  because  if  you  change  one  it  affects  the  other.  It  was  also  considered  that  the 
sum  of  the  variances  from  the  HL  and  BF  run  would  equal  the  variance  of  the  HLBF  run, 
which  is  not  the  case  either.  What  is  known  is  that  the  actual  uncertainty  due  to 
branching  fractions  must  be  between  BF  and  1  — ^hl_ 

°HLBF  °HLBF 

To  present  the  results  clearly,  a  count  of  the  number  of  isotopes  whose  ratio  is 
greater  than  a  set  amount  was  done.  As  an  example,  if  the  amount  was  0.5,  then 
1  — >0.5  or  gBF  >0.5  must  be  true  for  an  isotope  to  be  counted.  The 

°HLBF  °HLBF 

amount  0.5  would  count  all  isotopes  that  have  uncertainty  due  to  branching  fractions  as 
the  dominant  factor.  If  the  amount  was  0.2  then  all  isotopes  who  have  at  least  20%  of 
their  uncertainty  due  branching  fractions  is  counted.  This  count  of  isotopes  is  done  for 
both  ratios  gBF  and  1  — separately,  and  divided  by  the  total  number  of  isotopes 

°HLBF  °HLBF 

present  at  that  time.  This  gives  a  percentage  of  isotopes  that  meet  the  criteria.  As 
mentioned  in  the  previous  paragraph  the  actual  number  is  somewhere  between  the  two 
ratios.  To  take  this  into  account  the  two  counts  are  averaged  and  error  bars  show  the 
range  between  them,  Figure  8  shows  this  for  U238  high  where  uncertainty  due  to 
branching  fractions  dominate. 
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Figure  8.  Comparison  of  uncertainty  due  to  branching  fractions  ratios  as  a  percent  of  all  isotopes 

present  for  U238  high  and  <p  =  10%. 


Figure  8  shows  that  as  far  as  the  number  of  isotopes  goes,  the  difference  between 
the  two  ratios  is  typically  within  5%  of  each  other.  For  the  other  types  of  fission  the 
shape  and  magnitude  of  the  curves  are  similar  (±  2%).  The  one  with  the  largest 
difference  is  U238  fast,  whose  percent  of  isotopes  present  values  are  all  approximately 
5%  greater  at  all  times  shown. 

Figure  8  shows  the  point  at  which  either  uncertainty  from  half-lives  or  branching 
fractions  dominate.  To  get  a  better  idea  of  what  the  distribution  of  these  ratios  is  like  an 
empirical  CDF  was  done  at  an  early  time  and  late  time  and  is  shown  in  Figure  9. 
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Figure  9.  Comparison  of  empirical  CDFs  for  the  different  ratios  at  1  second  and  -1.5  years  for  U238 
high  and  <p  =  10%.  At  both  times  the  ratio  for  most  isotopes  is  either  very  close  to  zero  or  one.  This 
means  very  few  isotopes  have  uncertainties  from  both  half-lives  and  branching  fractions  that  are  of 

the  same  order  of  magnitude. 


The  empirical  CDFs  in  Figure  9  show  that  at  both  early  and  late  times  the  isotopes 
present  either  are  significantly  dominated  by  uncertainty  from  half-lives  or  uncertainty 
from  branching  fractions.  In  all  four  curves  the  ratio  values  go  from  less  than  0.1  to 
greater  than  0.9  in  less  than  0.2  on  the  CDF  axis.  This  means  less  than  20%  of  the 
isotopes  present  at  both  times  have  comparable  uncertainties  from  both  half-lives  and 
branching  fractions. 
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Another  way  of  showing  what  is  going  on  and  confirms  this  trend  over  all  times  is 


to  graph  multiple  curves  with  different  ratio  values,  this  is  done  in  Figure  10. 


Figure  10.  Comparison  of  multiple  uncertainty  ratio  values  as  a  percent  of  all  isotopes  present  for 

U238  high  and  <p  =  10%. 


Figure  10  shows  how  closely  spaced  the  ratio  values  are,  meaning  at  all  times  a 
majority  of  the  distribution  is  either  very  close  to  zero  or  one.  Again  this  supports  that 
for  most  isotopes  their  uncertainty  is  either  dominated  by  half-lives  or  branching 
fractions,  very  few  have  comparable  contributions  from  both.  Both  Figure  9  and  Figure 
10  look  at  only  U238  high  with  (p  =  10%.  Other  fission  types  and  tps  were  plotted  and 
showed  very  similar  results. 
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Branching  Fraction  Uncertainty  {(p)  Results 

Now  that  a  method  for  comparing  all  the  isotopes  present  at  a  given  time  has  been 
developed,  the  branching  fraction  uncertainty  ( cp )  can  be  varied  and  analyzed.  As  stated 
earlier  runs  were  done  with  (p  =  1%,  5%,  10%,  and  20%.  Ratio  values  of  0.5  were  chosen 
for  the  curves  because  it  shows  the  point  where  either  uncertainty  due  to  branching 
fractions  or  half-lives  dominate  the  overall  uncertainty.  Figure  1 1  shows  the  effect  of 
varying  <p. 


Figure  11.  Comparison  of  Branching  Fraction  Uncertainties  (<p)  with  ratio  values  of  0.5  as  a  percent 

of  all  isotopes  present  for  U238  high. 
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What  is  apparent  in  Figure  1 1  is  that  changing  <p  has  little  effect  on  changing  the 
percent  of  isotopes  it  dominates  at  a  given  time.  Going  from  cp  =  1%  to  (p  =  20% 
typically  only  increases  the  percent  of  isotopes  present  by  about  10%.  It  should  be  noted 
that  at  early  times  the  number  of  isotopes  present  is  around  800,  where  at  late  times  it 
dwindles  down  to  a  little  over  200.  The  real  driving  factor  is  time.  There  is  an  increase 
from  less  than  10%  at  less  than  one  minute  to  around  50%  of  isotopes  present  in  the 
month  timeframe.  The  fact  that  time  is  a  factor  makes  sense.  As  isotopes  decay  away 
and  become  very  long  lived  or  stable,  half-life  uncertainty  should  begin  to  decrease.  This 
is  because  half-life  uncertainty  is  only  a  factor  in  determining  when  an  isotope  will  be 
present.  Branching  fraction  uncertainty  is  a  factor  in  determining  what  isotope  will  be 
present.  So  as  time  passes  and  more  long  lived  or  stable  isotopes  are  present  the  effect  of 
half-life  uncertainty  decreases,  whereas  as  time  passes  and  more  opportunities  for 
isotopes  to  branch  occur  branching  fraction  uncertainty  should  increase.  What  was 
surprising  in  Figure  1 1  is  how  quickly  the  curves  rise  to  where  a  majority  of  the  isotopes 
present  are  dominated  by  branching  fraction  uncertainty. 

As  shown  in  Figure  9  from  the  previous  section  the  ratio  distributions  go  from 
zero  to  one  so  quickly.  This  means  changing  the  ratio  value  in  Figure  1 1  anywhere 
between  0.05  and  0.95  would  raise  or  lower  the  curves  by  only  a  few  percent.  This  also 
helps  explain  as  to  why  there  is  little  difference  when  cp  is  changed.  Because  most 
isotopes  are  either  significantly  dominated  by  half-life  or  branching  fraction  uncertainty 
changing  (p  would  have  little  effect.  The  driving  factor  is  not  how  big  <p  is,  but  the  fact 
that  it  is  present.  This  applies  when  analyzing  all  possible  isotopes  present  at  a  given 
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time.  If  individual  isotopes  are  analyzed  varying  cp  can  play  a  role  in  the  magnitude  of 
the  uncertainty.  The  ground  state  of  Snl29  is  shown  as  an  example  in  Figure  12. 


Figure  12.  Empirical  CDFs  of  the  Ground  State  of  Snl29  with  different  values  of  cp  from  U235 

thermal  at  10,000  seconds. 


As  <p  increases  the  empirical  CDF’s  range  increases,  showing  that  the  uncertainty 
goes  up,  while  the  mean  remains  constant.  The  standard  deviation  of  these  distributions 
is  shown  in  Table  3,  and  shows  that  as  <p  increases  the  standard  deviation  increases 
proportionally. 
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Table  3.  Quantity  standard  deviation  of  Snl29  ground  state  for  varying  <p  at  10,000  seconds  from 

U235  thermal  fission. 


9 

1% 

5% 

10% 

20% 

Standard  Deviation 

1.0192e-17 

5.0594e-17 

1.0319e-16 

2.0316e-16 

Increasing  cp  does  have  an  effect.  However,  for  most  isotopes  increasing  <p  as 
compared  to  the  contribution  of  uncertainty  from  half-lives  has  little  effect.  This  is 
because  the  uncertainty  from  half-lives  is  most  likely  either  much  greater  than  or  much 
less  than  the  uncertainty  from  branching  fractions. 

Limiting  Quantity  Fraction 

OAO 

Double  precision  in  Fortran  allows  numbers  to  get  as  small  as  10'"  .  However, 
this  is  much  smaller  than  what  is  needed  to  realistically  track  quantity  fractions  of 
isotopes.  If  there  are  approximately  1.45e“  fissions  per  kiloton  of  yield  Bridgman  [6], 
and  the  largest  potential  weapon  would  be  around  100  megatons,  this  means  a  quantity 
fraction  would  have  to  be  greater  than  6.9e'  to  be  present.  For  a  1  kiloton  weapon  this 
number  would  be  6.9e'“  .  Potential  collection  samples  would  have  to  be  much  larger  than 
this.  This  number  is  called  the  limiting  quantity  fraction  or  quantity  which  we  care  about 
or  can  measure.  It  is  then  used  as  a  cutoff  where  any  quantity  fraction  less  than  this 
limiting  quantity  fraction  is  not  counted  in  statistics.  The  limiting  quantity  fraction  can 
be  varied  to  see  how  it  affects  the  percentage  of  isotopes  present  that  are  dominated  by 
branching  fraction  uncertainties. 
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Figure  13.  Limiting  quantity  fraction  effect  on  percent  of  isotopes  dominated  by  branching  fraction 

uncertainty  for  U238  high  fission  and  tp  =  10%. 


Figure  13  shows  four  different  minimum  quantities  showing  the  line  for  percent 
of  isotopes  dominated  by  branching  fraction  uncertainty.  The  error  bars  are  omitted  in 
this  figure  so  the  reader  can  discern  the  small  differences  in  the  curves.  At  early  times 
there  is  very  little  difference  in  the  curves  and  all  are  well  within  each  other’s  error  bars. 
At  late  times  the  smallest  two  levels  are  still  right  on  top  of  each  other.  This  shows  there 
is  a  level  where  limiting  quantity  fraction  has  little  or  no  effect.  The  quantity  fraction  is 
so  small  that  most  isotopes  that  are  at  this  level  are  most  likely  rapidly  decaying  away  to 
zero.  The  interesting  thing  is  as  the  limiting  quantity  fraction  is  increased  the  curves 
grow  in  magnitude  as  time  increases.  The  purple  curve  representing  a  limiting  quantity 
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fraction  of  10"5  is  about  the  top  50%  isotopes  in  quantity.  Meaning  at  early  times  it  only 
shows  about  400  isotopes,  and  at  late  times  a  little  over  100  isotopes.  This  suggests  that 
isotopes  present  in  larger  amounts  are  more  likely  to  be  dominated  by  branching  fraction 
uncertainties  rather  than  half-life  uncertainties  at  later  times. 

Quantity  Fraction  Distributions 

One  of  the  program  outputs  is  each  isotope’s  quantity  fraction  mean  and  standard 
deviation  for  each  time  and  initial  quantity  value.  This  allows  for  comparison  of  isotopes 
across  times  or  different  initial  quantities  with  error  bars  that  represent  the  standard 
deviation.  In  the  absence  of  any  additional  information,  it  is  common  practice  to  assume 
that  these  statistics  are  normally  distributed.  The  goal  of  this  section  is  to  show  when  a 
normal  distribution  is  and  is  not  a  good  approximation  of  an  isotope’s  quantity  fraction. 

Analysis  of  all  isotopes  present  after  fission  for  all  times  would  be  overwhelming. 
Instead  individual  isotopes  that  displayed  a  particular  characteristic  across  multiple 
fission  types  or  timeframes  were  analyzed.  This  was  done  to  help  try  and  identify  trends 
or  attributes.  The  first  isotope  chosen  was  the  ground  state  of  Snl29.  It  was  chosen 
because  it  is  dominated  by  half-life  uncertainty  at  10,000  seconds.  Figure  14  shows 
empirical  CDFs  for  a  runs  with  branching  fraction  uncertainty  only,  half-life  uncertainty 
only,  and  with  both  a  branching  fraction  and  half-life  uncertainty.  Overlaid  with  dashed 
lines  are  normal  CDFs  with  means  and  standard  deviations  that  match  the  isotope  data. 
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Figure  14.  Empirical  CDFs  of  the  Ground  State  of  Snl29  from  U235  thermal  at  10,000  seconds  with 
<p  =  1%,  as  compared  to  a  Normal  CDFs  with  the  same  mean  and  standard  deviation. 


As  can  be  seen  in  Figure  14  when  half-life  uncertainty  is  considered  compared  to 
branching  fraction  uncertainty  it  increases  the  standard  deviation.  It  also  deviates  more 
from  a  normal  distribution  in  both  the  tails  and  around  the  peak.  This  carries  through  to 
the  run  when  both  uncertainties  are  considered.  Even  so,  all  three  runs  match  up  closely 
with  their  respective  normal  distributions.  The  last  thing  to  note  about  this  isotope  is  it 
just  starting  to  decay  away  to  zero,  and  within  one  more  decade  or  100,000  seconds  it 
will  be  gone. 
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Now  an  isotope  whose  uncertainty  is  dominated  by  branching  fractions  will  be 


examined.  Tel 24  is  very  long  lived  and  is  dominated  by  branching  fraction  uncertainty 
for  the  entire  timeframe  examined.  Figure  15  looks  at  the  distributions  at  early  time  0.1 
seconds  for  all  three  types  of  runs. 


0.9 


0.8 


0.7 


0.6 


Q  0.5 
O 


Branching  Fraction 
Uncertainty  Empirical  CDF 
Half-life  Uncertainty 
Empirical  CDF 

Half-life  and  Branching  Fraction 
Uncertainty  Empirical  CDF 
Branching  Fraction 
Uncertainty  Normal  CDF 
Half-life  Uncertainty 
Normal  CDF 

Haif-life  and  Branching  Fraction 
Uncertainty  Normal  CDF 


0.4 


0.3 


0.2 


0.1 


8.6  0.8 


1.2  1.4  1.6  1.8  2 

Quantity  Fraction 


2.2  2.4 


,.-io 

x  10 


Figure  15.  Empirical  CDFs  of  Tel24  from  Pu239  fast  at  0.1  seconds  with  <p  =  20%,  as  compared  to  a 
Normal  CDFs  with  the  same  mean  and  standard  deviation. 

Even  though  branching  fraction  uncertainty  is  the  dominant  factor  at  this  time,  the 
distributions  remain  very  close  to  normal  in  Figure  15.  Another  difference  between 
Figure  14  and  Figure  15  is  the  difference  in  <p.  In  Figure  14  <p  =  1%,  and  in  Figure  15  <p  = 
20%,  however  both  curves  with  branching  fraction  uncertainty  only,  match  up  closely 
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with  their  normal.  This  shows  (p  does  not  affect  the  shape  of  the  distribution.  There  was 
some  thought  that  since  branching  fractions  are  sampled  from  beta  distributions  when  the 
uncertainty  gets  large  it  could  deviate  more  from  normal  distributions  than  half-life 
uncertainty,  but  this  looks  like  it  is  not  the  case.  During  early  time  Tel24  is  going 
through  a  transient  and  is  increasing  in  quantity  as  seen  in  Figure  16  with  a  log  log 
quantity  versus  time  plot. 
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Figure  16.  Log  log  plot  of  the  quantity  of  Tel24  as  a  function  of  time  from  Pu239  fast  where  <p  = 
20%.  The  error  bars  represent  one  standard  deviation  from  the  mean. 
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The  transient  at  early  times  increase  the  quantity  by  several  orders  of  magnitude 
but  the  error  bars  stay  within  that  order  magnitude,  and  from  Figure  15  the  distributions 
remain  close  to  normal. 

The  next  isotope  analyzed  was  Xel39.  It  was  also  chosen  because  around  100 
seconds  its  quantity  is  very  stable  with  no  transients.  Figure  17  shows  CDFs  for 
branching  fraction  compared  to  half-life  and  branching  fraction  uncertainty. 
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Figure  17.  Empirical  CDFs  of  the  Xel39  from  U238  fast  after  100  seconds  with  <p  =  1%,  as  compared 
to  a  Normal  CDFs  with  the  same  mean  and  standard  deviation. 


Both  the  empirical  and  normal  CDFs  in  Figure  17  match  up  very  closely.  They 


are  the  best  fits  to  normal  distributions  thus  far.  When  half-life  uncertainty  is  added  in  it 
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can  be  seen  the  standard  deviation  and  uncertainty  increase  significantly.  This  means 
half-life  uncertainty  dominates.  The  half-life  uncertainty  alone  was  not  shown  in  Figure 
17  because  it  lies  directly  on  top  of  the  half-life  and  branching  fraction  uncertainty  CDFs. 
The  best  explanation  for  the  Xel39  CDFs  being  so  close  to  normal  is  that  it  is  not  going 
through  any  quantity  transients. 

Now  an  isotope  that  has  a  very  large  quantity  uncertainty  at  a  time  needs  to  be 
examined.  Figure  18  shows  the  quantity  of  Ni74  as  a  function  of  time  for  two  runs,  with 
only  half-life  uncertainty  as  well  as  half-life  and  branching  fraction  uncertainty. 
Branching  fraction  uncertainty  was  left  out  because  its  small  error  bars  cluttered  the  plot, 
and  added  nothing. 
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Figure  18.  Log  log  plot  of  the  quantity  of  Ni74  as  a  function  of  time  for  U235  fast  fission  with  <p  = 
5%.  The  error  bars  represent  one  standard  deviation  from  the  mean. 

As  Ni74  begins  rapidly  decaying  away  to  zero  the  two  data  sets  begin  to  diverge. 
This  shows  how  half-life  uncertainty  can  change  the  quantity  of  a  rapidly  decaying 
isotope  by  orders  of  magnitude.  This  is  also  shown  by  looking  at  the  error  bars  in  the  last 
three  time  steps  where  they  only  have  an  upper  bound.  This  also  corresponds  to  the  times 
when  its  standard  deviation  is  larger  than  the  mean.  Examining  this  further,  Figure  19 
takes  a  closer  look  at  the  empirical  CDFs  of  the  last  time  step  (46.4  seconds)  where  this 
isotope  is  present. 
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Figure  19.  Empirical  CDFs  of  Ni74  from  U235  fast  fission  at  46.4  seconds  with  <p  =  5%. 

Note  that  the  quantity  axis  in  Figure  19  is  logarithmic.  This  obviously  means 
these  distributions  are  not  normal,  but  more  closely  follow  some  kind  of  log  distribution. 
It  also  means  the  error  bars  at  these  late  times  are  not  accurate,  and  must  be  adjusted  to 
accurately  represent  this  logarithmic  distribution.  This  was  the  only  condition  found  that 
caused  the  quantity  distributions  to  differ  from  normal.  This  condition  is  when  half-life 
uncertainty  is  present  as  the  isotope  is  decaying  away  rapidly  to  zero.  The  empirical  CDF 
with  only  branching  fraction  uncertainty  was  also  examined  but  again  it  matched  a 
normal  distribution  as  in  the  two  previous  example  isotopes.  The  CDFs  from  an  earlier 
time  (1  second)  when  Ni74  is  not  rapidly  decaying  away  is  shown  in  Figure  20. 
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Figure  20.  Empirical  CDFs  of  Ni74  from  U235  fast  at  1  second  with  <p  =  5%. 


These  CDFs  match  up  much  closer  to  normal  distributions,  they  are  at  least  not 
logarithmic.  This  means  the  error  bars  in  Figure  18  at  early  times  are  accurate,  and 
shows  that  as  an  isotope  gets  close  to  decaying  away  its  quantity  distribution  moves  to  a 
logarithmic. 

It  also  shows  how  quantity  uncertainty  when  dominated  by  branching  fraction 
uncertainty  can  be  closely  approximated  by  a  normal  distribution.  This  was  not  always 
the  case  when  half-life  uncertainty  was  dominate,  especially  when  the  isotope  is  rapidly 
decaying  away. 
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V.  Conclusions 


The  objective  of  this  study  was  to  create  a  program  that  could  be  used  to  explore 
the  effects  of  branching  fraction  uncertainty  on  radioactive  decay  chains.  The  program 
leveraged  work  already  done  to  model  half-life  uncertainty  using  the  exponential 
moments  method  and  transmutation  matrices.  A  method  was  developed  for  random 
sampling  of  branching  fractions  using  a  Dirichlet  distribution  marginalized  by  beta 
distributions.  A  sampling  rejection  scheme  was  implemented  to  sample  these  beta 
distributions  parsimoniously  and  quickly. 

Because  branching  fraction  uncertainty  data  for  most  isotopes  is  unknown,  an 
unbiased  method  for  calculating  uncertainties  was  developed.  This  method  took  the 
branch  with  the  smallest  mean  of  each  isotope  and  made  its  relative  standard  deviation  a 
certain  percentage  (q>).  This  percentage  was  an  input  parameter  of  the  program,  and  was 
the  same  for  all  isotopes.  Nine  runs  with  1,000  Monte  Carlo  trials  each  were  completed. 
One  run  was  done  with  only  half-life  uncertainty  considered.  Two  runs  were  done  at 
with  (p  equal  to  1%,  5%,  10%,  and  20%  each,  one  being  with  only  branching  fraction 
uncertainty,  and  the  other  with  both  half-life  and  branching  fraction  uncertainty 
considered. 

The  standard  deviation  from  each  isotope  present  at  a  given  time  for  each  of  the 
different  runs  was  compared.  Using  ratios  of  these  uncertainties,  overall  uncertainty  from 
branching  fractions  was  compared  to  uncertainty  due  to  half-lives.  As  a  percentage  of 
isotopes  present  at  a  given  time  it  was  shown  that  uncertainty  due  to  half-lives  dominated 
most  isotopes  at  early  times  (<  1  month),  and  uncertainty  due  to  branching  fractions 
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dominated  at  times  greater  than  1  month.  This  makes  sense  because  half-life  determines 
more  of  when  an  isotope  is  present,  where  branching  fractions  determine  what  amount  of 
that  isotope  is  present.  As  time  passes,  uncertainty  due  to  half-life  decreases  while 
uncertainty  from  branching  fractions  increases  as  more  decays  occur  introducing  more 
opportunities  for  branching.  Increasing  cp  from  1%  to  20%  only  increased  the  number 
isotopes  dominated  by  branching  fraction  uncertainty  by  about  5%. 

It  was  also  shown  that  approximately  20%  of  isotopes  present  at  any  given  time 
have  uncertainties  of  the  same  order  of  magnitude  from  both  half-lives  and  branching 
fractions.  This  means  a  majority  of  isotopes  are  either  dominated  by  uncertainty  from 
half-lives  or  branching  fractions.  When  only  isotopes  present  in  greater  quantities  are 
considered,  the  percent  of  isotopes  present  dominated  by  branching  fractions  increases  at 
times  later  than  1  minute. 

The  quantity  distributions  about  the  means  can  be  considered  very  close  to  normal 
regardless  of  if  uncertainty  due  to  half-life  or  branching  fraction  dominates.  The  only 
known  exception  to  this  is  when  an  isotope  begins  rapidly  decaying  to  zero  on  a  log  time 
scale.  As  this  occurs  the  distribution  becomes  logarithmic. 

This  study  only  looked  for  general  guidelines  that  can  be  applied.  The  program  is 
flexible  and  the  data  outputs  sufficient  enough  to  study  particular  chains  or  isotopes  of 
interest  to  get  a  better  idea  of  branching  uncertainty  effects  on  them. 

Potential  Future  Work 

It  should  be  remembered  that  this  work  was  only  a  parameter  study,  with  all  the 
branching  fraction  uncertainties  being  applied  in  a  somewhat  equitable  manner.  The 
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most  logical  improvement  to  the  program  would  be  to  use  experimentally  determined 
branching  fraction  uncertainty  data.  If  enough  of  these  data  become  available  to 
incorporate,  it  should  be  remembered  that  if  the  variance  of  one  of  the  branches  for  an 
isotope  is  known,  the  variance  of  the  other  branches  is  already  determined.  This  means  if 
variances  for  multiple  branches  are  experimentally  determined,  a  method  for  weighting 
each  dataset  must  be  developed  so  as  to  give  an  accurate  Dirichlet  distribution.  Another 
improvement  would  be  to  include  uncertainty  of  the  initial  quantity  of  isotopes  from 
neutron  induced  fission.  If  this  is  done,  it  should  be  remembered  this  can  also  be 
represented  with  a  Dirichlet  distribution  because  it  is  categorical.  However  the  Dirichlet 
distribution  will  have  many  more  parameters,  upwards  of  821,  with  the  potential  for  very 
large  parameters  that  could  not  be  randomly  sampled  with  double  machine  precision. 

Improvements  to  the  efficiency  of  the  program  are  another  area  for  future  work. 
These  include  saving  out  the  all  the  decay  chains  to  an  array  so  they  would  not  have  to  be 
calculated  each  Monte  Carlo  run,  or  a  more  parsimonious  way  of  sampling  the  half-life 
uncertainty  distributions. 
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Appendix  A:  Design  Diagram  of  Decay  Program 
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Appendix  B:  Derivation  of  Standardized  Beta  pdf 


Probability  density  function  for  a  beta  distribution  with  parameters  a  and  /?. 

xa_1  (1  —  xY”1 
Betafett,W=  B(ap) 

where  B  (a, /?)  =  ta~1(l  — 

let  a  =  a  —  1,  and  b  =  ft  —  1 
with  a  domain  of  a  >  0,  b  >  0,  and  0  <  x  <  1 


(23) 


The  beta  distribution  pdf  is  now: 


Beta(x;  a,  b )  = 


xa(l  —  x)1 


B(a  + 1,6  +  1)  (24) 

To  standardize  the  pdf  the  max  value  must  be  known  (where  the  peak  is),  this  is  done  by 
finding  where  the  derivative  with  respect  to  x  is  equal  to  zero. 

(1  —  x)b~1xa~1(a(  1  —  x)  —  bx) 

Betafca'i0= - B(a  + 1,6  +  1) - =°  (25) 

The  only  way  Equation  (25)  can  be  equal  to  zero  is  when  (a(  1  —  x)  —  bx)  =  0.  Solve 
for  x: 


a  a  —  1 

X  a  +  b  a  +  —  2  (26) 

this  equation  is  equal  to  the  mode  of  a  Beta  distribution. 

Now  the  pdf  can  be  standardized  to  a  peak  value  of  one  by  dividing  by  the  peak 

value. 


Beta(x;  a,  b) 

Beta(^T6:a',,) 


((a  +  6)(  1  —  x) 

V  b 


(a  +  b)x 


\-  b)x\ 
a  ) 


(27) 
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This  is  a  fairly  simple  equation  that  is  not  only  standardized  to  a  peak  value  of  one,  but 
also  does  not  have  a  Beta  function  in  it.  It  can  be  further  simplified  to: 


Betas(x;  a,b ) 


(a  +  b)a+b 
aabb 


xa(l  —  x)b 


(28) 


The  issue  with  Equations  (27)  and  (28)  is  that  when  either  a,  b,  or  a  +  b  get  large  they  can 
numerically  overflow.  This  is  resolved  by  taking  the  log  of  each  side  and  leveraging 
logarithm  rules  so  Equation  (28)  becomes: 

Betacfx-a  b)  —  e(a+b^°S(a+b)-a\og(a)-b\og(b)+a\og(x)+b  log(l-x) 

>  >  s  ( OQ't 
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Appendix  C:  Derivation  of  Beta  pdf  Inflection  Points 


Take  either  Equation  (27)  or  (28)  and  find  its  second  derivative. 

/(a  +  b)(i  -  x)\  /(a  +  b)A  (a2(x  _  ty  +  {b_  l)bx2  +  a(x  _  1)(1  +  (2h  _  1)je)) 

Betas"(x;  a,b)  —  - - - - >  V  ’ - — — -  (30) 

( x  —  \)  X 

The  inflection  points  are  where  the  second  derivative  is  equal  to  zero.  The  only  way  for 
this  to  happen  in  Equation  (30)  is  for  the  third  factor  to  be  equal  to  zero. 

(a2(x  —  l)2  +  (b  —  1  )bx2  +  a(x  —  1)(1  +  (2  b  —  l)x))  =  0 

Solve  for  x: 


a2  +  a(b  —  1)  +  Jab(a  +  b  —  1) 


X  a2  +  (h  —  1  )b  +  a(2b  —  1)  (31) 

If  a  <  1  or  b  <  1,  then  one  or  both  inflection  points  will  be  outside  x’s  range  of  (0,1).  If 
this  is  the  case  then  the  rejection  scheme  only  uses  the  relative  slope  points  and  the 
constant  around  the  peak  of  the  standardized  pdf.  As  a  and  b  approach  one  in  the  minus 
case,  loss  of  precision  occurs.  This  is  taken  care  of  by  multiplying  Equation  (31)  by: 

/a2  +  a(b  —  1)  —  yjab{a  +  b  —  1)\  /a2  +  a(b  —  1)  +  A/ah(a  +  b  —  1)\ 

X  V  a2  +  (b  —  l)b  +  a(2b  -  1)  )  \a2  +  a(/,  _  i)  +  Jab(a  +  b-  1  )J 

If  the  numerator  and  denominator  are  expanded  and  terms  cancelled  Equation  (32)  is  the 
result. 


a2  —  a 

a2  +  a(b  —  1)  +  yjab{a  +  b  —  1)  (32) 

Equation  (32)  is  now  well-conditioned  and  can  be  used  for  the  minus  case  of 
Equation  (31). 
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Appendix  D:  Derivation  of  Inverting  CDF  for  a  line 


For  a  line  5  with  a  left  bound  of  xq  and  right  bound  of  jci,  let  s(.i'o)  =  sq,  s(xi)  =  ,V| 


then  the  equation  for  the  line  is: 


s(it)  =  s0(l  —  it)  +  stu 


where  u  — 


x  —  Xn 


X-,  -  Xn 


and  0  <  u  <  1 


Integrate  Equation  (33)  to  get  the  total  area  under  the  line. 


s(u)du  —■ 


Sn  +5! 


Then  divide  the  equation  for  the  line  by  Equation  (34)  and  integrate  to  get  its  CDF. 


"  2  s(u) 

- du  — 

J  s0  +  s1 


2s0u  —  s0u2  +  sxul 
Sn  +Si 


Set  Equation  (35)  equal  to  <f,  a  random  number  between  [0,1)  and  solve  for  u. 


u(0  = 


So  ±VS0  ~S0?  +  si( 


Sn  Si 


To  find  the  correct  root  solve  u(  1),  and  u( 0)  for  each  case.  For  the  plus  case  (0)  = 
—  ,u(  1)  =  S°+Sl  ,  for  the  minus  case  u(0)  =  0,  u(  1)  =  1.  The  minus  case  is  correct 

So -Sl  So -Si 

for  the  CDF.  Expanding  Equation  (36)  and  factoring  to  remove  the  s0  —  from  the 
denominator  results  in: 


u(0  = 


(sQ  +  gi)£ 
s0  +  a/so(1~0 
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The  (1  —  <f)  in  the  denominator  will  not  result  in  loss  of  precision  because  <f  lives  in 
[0,1).  Now  u  can  be  randomly  sampled  from  the  CDF  of  the  line  s.  Because  of  the 
method  developed  with  four  linears  and  a  constant  two  of  the  sampling  function  lines  will 

have  either  so  =  0,  or  si  =  0.  When  so  =  0,  Equation  (37)  simplifies  to  y[%,  and  when 
si  =  0  it  simplifies  to  ^/(l  +  ^  1  —  £) . 
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Appendix  E:  Derivation  of  Standardized  Beta  pdf  Relative  Slope 


Start  with  either  Equation  (27)  or  (28)  from  Appendix  B,  and  taking  its  derivative 
you  get  an  equation  for  the  slope  of  the  standardized  beta  pdf: 


Betas'(x;  a,b ) 


((a  +  m  -  0)t((o (a(1  _  x)  _  fa) 


(1  —  x)x 

Then  divide  the  slope  by  the  standardized  pdf  to  get  the  relative  slope. 

Betas'(x;  a,  b )  a(x  —  1)  +  bx  a  b 

=  SR  —  7  —  —  —  + 1 


Beta s(x;a,  h)  (x  —  l)x  x  ’  (x  —  1) 

where  Sr  is  the  relative  slope.  Solving  for  x  results  in  two  solutions: 

a  +  b  +  SR  ±  yj(a  +  b  +  SR)2  —  4 aSR 


x  - 


2SK 


(38) 


(39) 


(40) 


of  which  the  minus  is  correct  because  it  yields  solutions  within  the  domain  0  <  x  <  1. 
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Appendix  F :  Prioritized  List  of  Model  and  Code  Improvements 

1.  Parallelization:  Because  the  program  relies  on  Monte  Carlo  methods, 
there  is  great  potential  that  parallelization  could  reduce  runtimes.  However  the 
auto  parallelization  feature  in  the  Intel  compiler  actually  increased  runtimes.  An 
analysis  of  the  code  by  someone  experienced  in  parallelization  techniques  is 
needed.  With  today’s  multi-processor  computers  this  could  lead  to  the  greatest 
reduction  in  runtime. 

2.  Incorporate  known  branching  fraction  uncertainties:  Some  branching 
fraction  uncertainties  are  documented,  and  more  may  be  done  in  the  future.  A 
method  that  incorporates  these  uncertainties  while  still  adhering  to  the  model  laid 
out  in  this  study  is  needed. 

3.  Calculate  decay  chains  only  once:  Currently  decay  chains  for  all  the 
isotopes  are  calculated  for  each  time  of  each  Monte  Carlo  trial.  These  decay 
chains  do  not  change,  making  this  method  redundant.  The  large  number  of  decay 
chains  may  make  it  infeasible  to  save  to  an  array,  or  too  inefficient  to  write  to  a 
file  and  read  back  in  every  time.  However,  these  or  other  methods  need  to  be 
explored  for  efficiency  sake. 

4.  Improve  sampling  of  half-life  uncertainties:  The  method  used  is  not 
parsimonious  in  random  numbers,  and  is  not  realistic  for  a  Gaussian  distribution 
(when  an  isotope  is  quoted  as  having  a  +a  and  -a).  It  was  recommended  by  NIST 
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to  splice  the  two  halves  of  the  different  Gaussians  together.  If  this  is  done  you  get 


a  distribution  similar  to  the  one  in  Figure  21. 


Figure  21.  pdf  of  two  Gaussians  with  different  standard  deviations  spliced  together.  The  mean  is  0.5. 


Approximately  200  isotopes  in  the  current  NuDat  data  have  two  different  a 
values.  A  more  appropriate  distribution  is  needed  that  does  not  have  a  step  value 
at  the  mean,  and  satisfies  both  a  values. 
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5.  Incorporate  uncertainty  in  initial  quantities:  This  is  another  piece  in  the 
physics  of  radioactive  decay  that  needs  to  be  explored  as  a  significant  source  of 
uncertainty.  As  an  example  fission  fragment  yields  range  across  many  orders  of 
magnitude  (10  '  to  10"  ).  Because  of  this,  uncertainty  in  these  yields  could  have 
an  impact  in  decay  isotope  quantities. 

6.  Account  for  spontaneous  fission:  For  this  study  the  fraction  of  an  isotope 
that  can  spontaneously  fission  was  ignored.  This  could  be  done  because  no 
isotopes  that  result  from  neutron  induced  fission  can  spontaneously  fission.  This 
does  need  to  be  accounted  for  to  make  the  code  complete. 

7.  Improve  the  Beta  distribution  rejection  scheme:  The  current  rejection 
scheme  has  an  average  acceptance  rate  of  around  63%.  If  the  goal  of  the  user  is 
run  more  Monte  Carlo  trials  and  only  a  few  times,  then  improving  this  scheme 
would  result  in  a  noticeable  reduction  in  runtime. 

8.  Transmutation  matrix  inversion:  Inverting  the  T-matrices  would  allow  the 
user  to  calculate  quantities  present  at  earlier  times  with  uncertainty.  This  is  useful 
in  nuclear  forensics,  where  the  quantity  of  an  isotope  is  measured  at  a  time  after  a 
nuclear  event.  The  need  is  to  be  able  to  plug  this  amount  into  a  program  that  will 
give  the  isotopes  initially  present  with  uncertainty. 
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